Scientific journal
Fundamental research
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,674

SIMULATION OF NON-LINEAR DYNAMICAL SYSTEMS WITH PARALLEL METHODS OF NUMERICAL INTEGRATION

Butusov D.N. 1 Ostrovskiy V.Y. 1 Krasilnikov A.V. 1
1 Saint-Petersburg State Electrotechnical University
The paper discusses the new class of ordinary differential equations solvers, based on parallel methods of numerical integration, known as D-methods. The application problems of such methods are considered for a case of non-linear dynamical systems computer simulation. The algorithms of parallel ODE integration methods and the series of computational experiments in NI LabVIEW development system are represented. The test examples of non-linear flat pendulum and Lorenz «strange attractor» system are given, solved by new parallel numerical integration methods. An evaluation of computational efficiency and stability of the obtained parallel method modifications is provided, compared to the classical Runge-Kutta methods.
non-linear dynamical system
numeric integration
ODE Solver
parallel computing
dynamical system simulation
1. Butusov D.N. Sintez i issledovanie apparatno-orientirovannyh chislennyh metodov integrirovanija v srede LABVIEW // Materialy XIII Mezhdunarodnoj nauchno-prakticheskoj konferencii «Inzhenernye i nauchnye prilozhenija na baze tehnologij National Instruments» NI Days 2014, M.: MTUSI, 2014.
2. Butusov D.N., Karimov A.I., Karimov T.I., Dolgushin G.K. Semejstvo apparatno-orientirovannyh metodov chislennogo integrirovanija // Sovremennye problemy nauki i obrazovanija. 2014. no. 4.
3. Zhukov K.G. Algoritm realizacii parallel’nyh vychislenij po formulam chislennogo integrirovanija Runge-Kutta // Nauchno-tehnicheskie vedomosti sankt-peterburgskogo gosudarstvennogo politehnicheskogo universiteta. Informatika. Telekommunikacii. Upravlenie. SPBGPU, 2011. T. 6.2, no. 138. pp. 143–149.
4. Fel’dman L.P. Parallel’nye algoritmy modelirovanija dinamicheskih sistem, opisyvaemyh obyknovennymi differencial’nymi uravnenijami // Jelektronnoe modelirovanie. 2004. t. 26, no. 1. pp. 19–30.
5. Hajrer Je., Nersett S., Vanner G. Reshenie obyknovennyh differencial’nyh uravnenij. Nezhestkie zadachi. – M.: Mir, 1990.

По мере развития средств вычислительной техники становится очевидной тенденция к широкому применению параллельных вычислений. Одной из основных задач, решаемых на ЭВМ, остается численное моделирование динамических систем, описываемых системами дифференциальных уравнений. Особый интерес представляет машинное решение нелинейных дифференциальных уравнений, как правило, не поддающихся решению аналитическим способом.

В большинстве инструментальных пакетов моделирования основным математическим аппаратом численного решения ОДУ остается семейство методов Рунге – Кутты. Эти хорошо изученные и зарекомендовавшие себя алгоритмы тем не менее имеют ряд существенных недостатков с точки зрения дальнейших перспектив развития моделирующих систем. В частности, алгоритмы методов Рунге – Кутты представляют собой жестко рекуррентные многостадийные итерационные процессы, плохо поддающиеся распараллеливанию. В последнее время возрос интерес к параллельным численным методам интегрирования, что отражено в большом количестве публикаций на данную тему, например [2, 3, 4]. При этом практически все известные авторам работы рассматривают варианты распараллеливания явных методов Рунге – Кутты или же непосредственно самих моделей динамических систем. При таком подходе сохраняются все недостатки явных методов Рунге – Кутты, связанные с их малой численной устойчивостью, накладывающей серьезные ограничения при моделировании жестких и разноскоростных систем. В настоящей работе рассматривается пример использования численных методов 2 и 4 порядка алгебраической точности, являющихся параллельными модификациями метода Эйлера и не относящихся к группе методов Рунге – Кутты, для численного моделирования нелинейных динамических систем.

Плоский нелинейный маятник

Классическая задача математического моделирования поведения плоского нелинейного маятника описывается следующим дифференциальным уравнением:

butusov01.wmf (1)

где l – длина нити между шаром и подвесом; m – масса шара; θ – угол между текущей и устойчивой позицией шара; τ – внешняя сила, воздействующая на шар.

Произведем замену переменных состояния θ и butusov02.wmf на z1 и z2 соответственно.

В нормальной форме Коши уравнение (1) примет вид

butusov03.wmf (2)

где butusov04.wmf butusov05.wmf

Алгоритм классического метода Рунге – Кутты второго порядка алгебраической точности (РК2) для решения системы (2) может быть описан следующей последовательностью действий:

butusov06.wmf

butusov07.wmf

butusov08.wmf

butusov09.wmf

butusov10.wmf

butusov11.wmf

butusov12.wmf

butusov13.wmf (3)

где h – величина шага интегрирования.

Классический алгоритм содержит 18 операций, необходимых для расчета решения на каждом шаге интегрирования, из них 8 умножений, 8 сложений и 2 операции взятия синуса.

Реализуем алгоритм (3) как виртуальный прибор среды моделирования NI LabVIEW. Лицевая панель – интерфейс виртуального прибора – показана на рис. 1.

pic_7.tif

Рис. 1. Ввод параметров моделирования плоского нелинейного маятника на лицевой панели виртуального прибора

Рассмотрим поведение переменных состояния и фазовый портрет дискретной системы (3), смоделированной численным методом РК2. На рис. 2 представлен график зависимости расчетного угла отклонения маятника от времени, а также фазовый портрет системы при шаге интегрирования 0,1 c и времени моделирования 45 с.

Из рис. 1 видно, что накапливающаяся ошибка решателя приводит со временем к неустойчивости дискретной системы и расходящемуся решению, что является следствием малой численной устойчивости явных методов Рунге – Кутты.

Рассмотрим аналогичный решатель для предлагаемого параллельного метода второго порядка точности [2], называемого методом второго порядка с диагональной коррекцией (далее Д2). Алгоритм метода при решении системы (2) имеет вид

butusov14.wmf

butusov15.wmf

butusov16.wmf

butusov17.wmf

butusov18.wmf;

butusov19.wmf, (4)

где h – величина шага интегрирования.

pic_8.tif

Рис. 2. Поведение системы (2) во временной области и фазовом пространстве при моделировании методом Рунге – Кутты 2

При этом можно заметить, что расчет z11 и z21 может быть выполнен независимо от расчета z12 и z22. При этом образуются два параллельных интегрирующих процессора. Усреднение значений переменных состояния, полученных на предыдущем шаге и являющихся выходными значениями на текущем, выполняются в составе параллельной части решателя (рис. 3) в том случае, если не требуется подача усредненного значения.

pic_9.tif

Рис. 3. Сравнение блок-диаграммы решателя системы (2) на основе метода РК2 (слева) и двухпроцессорного решателя на основе параллельного метода D2

Таким образом, алгоритм метода Д2 содержит 18 операций на шаге интегрирования, из которых 8 умножений, 8 сложений и 2 операции взятия синуса. Если при этом учесть параллельную структуру решателя, то количество последовательно выполняемых действий сокращается до 7 операций: из них 3 умножения, 3 сложения и 1 операция взятия синуса. Расчетное ускорение вычислений по сравнению с методом РК2 составляет около двух раз, причем можно предположить, что ускорение будет увеличиваться с ростом порядка моделируемой системы дифференциальных уравнений.

Кроме быстродействия и алгебраической точности при компьютерном моделировании динамических систем важную роль играет такая характеристика, как численная устойчивость метода. Известные аналитические способы определения устойчивости численных методов интегрирования предназначены для оценки устойчивости РК-методов, которые могут быть записаны в виде таблицы Бутчера. Однако Д-методы невозможно представить в виде таблицы Бутчера, т.к. они не опираются на обобщенную формулу Рунге – Кутты [5]. Это требует разработки другого подхода к оценке устойчивости Д-методов. Оценим устойчивость численного метода Д2 экспериментально, сравнив характеристики полученной с его помощью дискретной модели с характеристиками модели, построенной по методу РК2.

На рис. 4 представлен график поведения переменной состояния системы (2) и фазовый портрет при тех же параметрах моделирования (рис. 1), но при решении методом Д2.

pic_10.tif

Рис. 4. Поведение системы (2) во временной области и фазовом пространстве при моделировании параллельным методом Д2

При сравнении рисунков (2) и (4) становится очевидным, что модель, построенная по методу Д2, остается устойчивой, на всем интервале времени моделирования сохраняя заданную погрешность расчета значений переменных состояния, в отличие от метода Рунге ‒ Кутты второго порядка. Напомним, что при этом метод Д2 имеет в два раза меньше последовательно выполняемых операций на шаге.

Задача моделирования аттрактора Лоренца

Интересным направлением в моделировании нелинейных динамических систем являются системы с управляемым хаосом, ярким примером которых является известная система дифференциальных уравнений Лоренца, называемая в некоторых источниках [5] также «странным аттрактором Лоренца» (5).

butusov20.wmf (5)

Моделирование подобных систем предъявляет дополнительные требования к используемым численным методам, делая малопригодными методы второго порядка вследствие их малой области устойчивости и низкой точности решения. Поэтому для моделирования системы (5) используем параллельный метод четвертого порядка алгебраической точности Д4, описанный в работе [2]. Фазовый портрет системы (5) в трех измерениях, полученный путем моделирования со следующими параметрами: tмод = 80 с, σ = 10, r = 28, b = 2,6667, h = 0,01 с, представлен на рис. 5.

pic_11.tif

Рис. 5. Трехмерный фазовый портрет системы (5)

«Странный аттрактор» Лоренца, моделирование численным методом Д4.

При значении параметра r ≤ 13 значения переменных состояния системы сходятся к установившимся значениям, которые можно найти аналитически:

butusov21.wmf (6)

Используем данное свойство системы для оценки погрешности ее моделирования предлагаемым численным методом. Максимальная абсолютная погрешность моделирования системы (5) методом Д4 относительно точных значений установившегося решения при параметрах моделирования: tмод = 40 с, σ = 10, r = 13, b = 2,6667, h = 0,01 с составляет Eуст = 1,3506∙10–7, что сопоставимо с погрешностью явного метода Рунге – Кутты 4 порядка (Eуст = 1,3509∙10–7). Это подтверждает пригодность предложенного численного метода для численного моделирования динамических систем с управляемым хаосом. При этом, обладая свойством вычислительного параллелизма, метод Д4 обеспечивает как минимум двукратное ускорение процесса моделирования по сравнению с методом РК4 при реализации на ЭВМ с параллельной архитектурой [1].

Заключение

В статье показано, что предложенный параллельный численный метод интегрирования второго порядка алгебраической точности имеет бо́льшую численную устойчивость при моделировании нелинейных систем, находящихся на грани устойчивости, по сравнению с классическим методом Рунге – Кутты 2. При реализации решателя по методу Д2 на вычислителях с параллельной архитектурой также достигается двукратное ускорение моделирования за счет сокращения числа последовательных операций на шаге интегрирования. Приведен пример моделирования нелинейной хаотической динамической системы параллельным численным методом 4 порядка алгебраической точности, проведена оценка погрешности полученного решения в установившемся режиме. В качестве дальнейшего развития изложенных в статье идей предполагается синтез Д-методов более высоких порядков алгебраической точности, а также адаптация подобных методов к работе с переменным шагом интегрирования. Результаты исследования могут быть использованы при моделировании нелинейных динамических систем на ЭВМ с параллельной архитектурой, в том числе в реальном масштабе времени, а также в составе пакетов компьютерного моделирования.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований в рамках договора № НК 14-01-31277\14 от 06.02.2014 г.

Рецензенты:

Авдеев Б.Я., д.т.н., профессор кафедры информационно-измерительных систем и технологий, ФГАОУ ВО «Санкт-Петербургский государственный электротехнический университет «ЛЭТИ» им. В.И. Ульянова (Ленина)», г. Санкт-Петербург;

Анисимов В.И., д.т.н., профессор кафедры систем автоматизированного проектирования, ФГАОУ ВО «Санкт-Петербургский государственный электротехнический университет «ЛЭТИ» им. В.И. Ульянова (Ленина)», г. Санкт-Петербург.

Работа поступила в редакцию 28.12.2014.