CC BY
134
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 • 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Сравнение современных решателей жестких систем
обыкновенных дифференциальных уравнений с решателями Си
библиотеки SADEL
# 08, август 2012
DOI: 10.7463/0812.0445558
Жук Д. М., Маничев В. Б., Сахаров М. К.
УДК.519.622.2
Россия, МГТУ им. Н.Э. Баумана zhuk@bmstu.ru manichev@bmstu.ru max. sfn90@gmail. com
Введение
В математическом ядре любой CAE (Computer Aided Engineering) системы реализованы различные методы численного интегрирования обыкновенных дифференциальных уравнений (ОДУ) и дифференциально-алгебраических уравнений (ДАУ). Однако, большинство инженеров-проектировщиков не являются специалистами в области численных методов и зачастую не могут оценить адекватность результатов, выдаваемых системой. Поэтому, при компьютерном моделировании и анализе на основе дифференциальных уравнений высокие требования предъявляются к эффективности применяемых численных методов. Эффективность численных методов включает в себя два аспекта:
• достоверность и точность получаемых решений;
• затраты машинного времени и оперативной памяти на получение решений.
В данной работе авторы акцентируют внимание на проблеме получения достоверного и точного решения систем ОДУ-ДАУ. Известно, что решение жестких систем ОДУ традиционно вызывает затруднения у большинства
численных методов. В наших работах [1-3] было показано, что для гарантии получения корректного и достоверного решения жестких систем ОДУ численный метод их решения должен быть AL-устойчивым, т.е. абсолютно устойчивыми строго в левой (Absolute Left) комплексной полуплоскости. AL-устойчивые методы 2-го и 4-го порядков точности были реализованы в программе DMAN [4].
Главным недостатком реализации вычислительных методов в известных программных комплексах, является возможная выдача ошибочных результатов для жестких систем ОДУ без предупреждения пользователей об их недостоверности. В рамках этой работы жесткой будем считать систему ОДУ со степенью жесткости более 106.
Работа посвящена тестированию различных реализаций методов численного интегрирования, предназначенных для решения жестких систем ОДУ. Предложенный авторами набор тестовых задач был использован для тестирования реализаций методов численного интегрирования в системах MATLAB, Maple, Mathematica, а также в библиотеке NAG и последующего сравнения результатов с решателями библиотеки SADEL (Sets of Algebraic and Differential Equations solvers Library), разработанной авторами [1].
1. Алгоритмы решения систем ДАУ в библиотеке SADEL
В существующих программах-решателях систем ОДУ-ДАУ используются три классические постановки задачи.
1) Нормальная форма Коши
(1)
Где X - вектор дифференцируемых переменных системы размерностью т, F -
вектор-функция правых частей, t - независимая переменная (обычно - время). В
координатный базис входят только переменные состояния. Такая постановка
задачи называется явной. Также заданы начальные условия Х0 = А'(0)и отрезок
интегрирования t = [0, Тк], где заданное значение.
2) Дифференциально-алгебраическая форма разных индексов
VG(XtY) = О,
Где У - вектор алгебраических переменных размерностью к, С - вектор-функция размерностью к. В данном случае в базис входят как дифференцируемые переменные, так и алгебраические. Такой базис называется полным координатным базисом, а постановка задачи - полуявной. Также заданы начальные условия Х0 = АЛ(0), согласованные начальные условия У0 = У(0) и отрезок интегрирования С = [0,7^]. Начальные условия алгебраических переменных являются согласованными в начальный момент времени, если они удовлетворяют алгебраическим уравнениям.
Отметим, что индекс дифференцирования системы есть наименьшее число аналитических дифференцирований, требующихся для того, чтобы из уравнений (2) получить систему (1). При этом каждое дифференцирование понижает индекс на единицу.
3) Дифференциально-алгебраическая форма в полном координатном базисе
Где G - вектор-функция размерностью т + к. Заданы начальные условия Х0 = Х(0), согласованные начальные условия У0 = У(0) и отрезок интегрирования t = [0, rfc]. Такая постановка задачи называется неявной.
Задача (3) была поставлена в полном координатном базисе Л. Петзольд в 1982 г. Это послужило основой для разработки знаменитой программы DASS (Differential Algebraic Systems Solver) на основе метода BDF (Backward Differential Formula).
Задача решения систем ДАУ в постановке (3), но в расширенном координатном базисе была поставлена в статье [5]. В данном случае в базис помимо дифференцируемых и алгебраических переменных включаются сами
(3)
производные дифференцируемых переменных. Этот базис используется в решателе систем ДАУ в библиотеке БЛОБЬ.
Для решения жестких систем ОДУ были разработаны формулы Б-стадийных ^¿-устойчивых неявных методов интегрирования систем ДАУ вида (3) в расширенном координатном базисе [6]. Формулы имеют вид
/с1Х \
АХ АХ
Ху1 Хг9 ! £ 1, ... , 5",
п
^Л
(4)
где 5 - число стадий, кп - м-й шаг интегрирования, с^ параметры
метода.
На основе этих формул были разработаны и реализованы три метода интегрирования [4]:
• М1 - А-устойчивый неявный метод Эйлера первого порядка точности;
• М2 - АЬ-устойчивый неявный метод второго порядка точности;
• М3 - АЬ-устойчивый неявный метод четвертого порядка точности.
На рисунках 1, 2 проиллюстрирована работа алгоритмов, реализующих эти методы. Методы основаны на совместном решении систем нелинейных
алгебраических уравнений С (х, У, ^^ = 0 и систем линейных алгебраических
уравнений Н = О? сформированных на соответствующих стадиях работы
алгоритмов [7].
Рис. 1. Схема алгоритмов неявных одно стадийных одношаговых методов Эйлера и трапеций (М1 и М2 - число стадий 5 = 1)
Рис. 2. Схема алгоритма неявного двух стадийного одношагового метода (М3-
число стадий 5 = 2)
В библиотеке БДОБЬ реализовано три метода решения систем линейных алгебраических уравнений (СЛАУ) [2]. Применение этих методов гарантирует решение с высокой точностью (15 верных значащих цифр) даже плохо обусловленных систем, с числом обусловленности многим больше единицы. Верное решение СЛАУ на каждой итерации алгоритмов интегрирования систем ОДУ-ДАУ способствует верному решение исходной задачи [8].
2. Алгоритмы решения жестких систем ОДУ в современных
решателях
Для решения жестких систем ОДУ традиционно используются неявные методы, так как они не имеют ограничения на величину шага интегрирования снизу, при условии, что в программной реализации формул неявного метода нет арифметических операций деления на шаг [9]. Одним из таких методов является метод БОБ.
Фактически, метод БОБ - это семейство неявных многошаговых методов численного интегрирования, применяемых для решения жестких систем ОДУ или систем ДАУ [10]. В качестве примера рассмотрим задачу (1). Общая формула для методов БОБ может быть записана в виде
(6)
где /г- это шаг интегрирования, ак_^ - коэффициент. Величина^ определяет порядок полинома аппроксимации и число интерполируемых точек до 1п+1. К примеру,
4
3
к — 2: Xn+i ~ ^Х-п
1 2
X-
71+1
18 9 2 6
11 ~ ll^71"1 + Vi^11'2 + ^ 11 ^^n*1^71-*-1^'
(7)
(8)
(9)
При к = 1 получаем неявный метод Эйлера, при к = 2 - метод BDF-2 и т.д. Отметим, что с ростом к область устойчивости метода уменьшается.
Метод BDF реализован практически во всех современных математических средах и библиотеках. В данной работе были протестированы реализации метода в ведущих программных пакетах-Wolfram Mathematica и MATLAB Simulink. Помимо этого была протестирована реализация метода в подключаемой библиотеке NAG, написанной на языке Си.
В программном пакете Maple отсутствует реализация метода BDF. Вместо нее для решения жестких систем ОДУ-ДАУ рекомендовано использование метода Розенброка (Rosenbrock) [11]. Метод Розенброка является L-устойчивым, он обеспечивает хорошую относительную точность Т = 10~4 (Tolerance) при небольшой размерности решаемой системы (iV < 10). Рассмотрим задачу (6). Метод Розенброка ищет решение системы в виде
(10)
где поправка Кь находится из решения системы 5 линейных уравнений. Система имеет вид
(11)
cLF
где--матрица Якоби. Коэффициенты у, с^а^/уц - константы не зависящие
dt ' ^ . i>~ij>rij
от задачи. Если у = уц = 0, то получаем метод Рунге-Кутта (Runge-Kutta).
3. Набор тестовых задач ОДУ-ДАУ
Набор тестовых задач состоит из четырех систем ОДУ. Эти системы имеют известное асимптотическое или аналитическое решение. Ниже приведено краткое описание данных тестов.
Тест 1, 1а - это классическая задача Ван дер Поля[12]. Тест представляет собой жесткую систему ОДУ второго порядка, где ¡х - параметр жесткости:
(12)
Сравнение программ было проведено для часто встречающихся на практике значений параметра жесткости/: = 106 (Тест 1) и /I = 109 (Тест 1а).
Тест 2, 2а - это математическая модель сверхколебательной электронной схемы. Модель представляет собой нежесткую система ОДУ пятого порядка с многопериодным аналитическим решением - тест Маничева [9]:
(13)
Д, = Ко = к
г>
£ = к
и'
к г
^ 1
к, =
к-^ * к.
и
"и
к = — г к,-
^ (Тест 2), а также для значений (Тест 2а), которые вызывают
Сравнение решателей было проведено для часто встречающихся на практике параметров масштабных коэффициентов по времени, току и напряжению - к{ = 1, к^ = 1 ,ки = 10 параметров к{ = 1СГ104,А:; = 1, ки = 1 затруднения у большинства современных решателей систем ОДУ.
Тест 3 представляет собой нелинейную жесткую систему ОДУ второго порядка, имеющую локально неустойчивое решение и составленную по примеру теста Ван дер Поля - тест Скворцова [13]:
Сравнение решателей было проведено для часто встречающегося на практике значения параметра жесткости ц = 106,
Тест 4 представляет собой нелинейную жесткую система ОДУ второго порядка для моделирования процессов работы лазера - тест Евстифеева [14]:
(15)
г = 0,016, Значения параметров данной системы соответствуют параметрам используемого на практике лазера.
4. Результаты сравнительного тестирования решателей
ОДУ-ДАУ
Тестирование проводилось для параметров решателей, установленных по умолчанию. Все решатели сравнивались при значении относительной точности 7 = Ю-", В таблице 1 приведены результаты сравнения решателей на предложенном наборе тестовых задач.
Таблица 1. Сводная таблица результатов тестирования решателей ОДУ-ДАУ
Тесты\Решатели ОДУ SADEL C-Library 2010 Методы M2,M3 Maple 2007 МетодРунг е-Кутта-Розенброка MATLAB 2007 Метод BDF NAG C-Library 2012 Метод BDF Mathematica 8 Метод BDF
Тест 1. Уравнения Ван дер Поля MU=106 + - + + -
Тест 1а. Уравнения Ван дер Поля MU=109 + - - - -
Тест 2. Высокодобротный фильтр Ы=1, Ы=1, кы=0.01 + + - - +
Тест 2а. Высокодобротный фильтр Ы=10-104, кы=1, кл=1 + - - - +
Тест 3. Локально неустойчивая система ОДУ MU=106 + - - - +
Тест 4. Моделирование свечения лазера + - + + +
Знак минус означает, что при стандартных настройках данных решателей авторам не удалось получить верное решение. В подавляющем большинстве случаев неверное решение не сопровождалось диагностическим сообщением о возможной ошибке. Из таблицы видно, что лишь решатели библиотеки SADEL успешно справились со всеми тестовыми системами.
Необходимо отметить, что реализация метода BDF в Wolfram Mathematica и реализации в MATLAB Simulink и NAG C-Library справились с разными тестами. Тот факт, что реализация в Mathematica справилась с большим количеством тестов, свидетельствует о модификации метода BDF разработчиками Wolfram.
При дополнительных настройках параметров решателей удалось существенно улучшить результаты тестирования. Например, при искусственном увеличение точности решатель библиотеки NAG не справился лишь с тестом 2а, а техника комбинирования явного метода Рунге-Кутта и неявного метода Эйлера при повышенной точности в пакете Wolfram Mathematica позволила успешно решить весь набор тестовых задач. Однако такие настройки требуют от
инженеров-проектировщиков широких знаний в области методов численного интегрирования систем ОДУ и хорошего владения используемыми математическими пакетами. Так как серьезной проблемой является почти полное отсутствие диагностических сообщений о возможных ошибках, то инженер должен также уметь оценить адекватность полученного решения. Решатели библиотеки SADEL позволяют решать любые системы ОДУ-ДАУ, выдавая гарантировано верный результат в широком диапазоне значений параметров жесткости, что существенно облегчает жизнь простому инженеру.
На рисунках 3-6 представлены некоторые результаты сравнения решателей ОДУ-ДАУ.
Рис. 3. Решение теста 2а полученное с помощью библиотеки NAG при Т = 10
-з
Рис
12 ООО
1.5
0.5
О
О 0.5 1 13 2 2.5 3
t
Рис. 5. Решение теста 3 полученное с помощью библиотеки NAG при Т = 10~3
О 0.5 1 1.5 2 2.5 3
t
Рис. 6. Правильное решение теста 3
Заключение
Проведенное сравнение современных решателей жестких систем ОДУ-ДАУ с решателями разработанной авторами библиотеки SADEL позволило сделать следующие выводы.
1) Даже лучшие CAE подсистемы на сегодняшний день не способны гарантировать качественно верное решение различных систем ОДУ. Также проблемой является отсутствие диагностических сообщений в случае некорректного решения.
2) Если программа-решатель систем ОДУ-ДАУ не решит достоверно и точно хотя бы одну из приведенных в статье тестовых задач, то этот решатель не следует использовать для моделирования динамических систем, математическая модель которых представляет собой жесткую систему ОДУ-ДАУ.
Новые научные результаты предполагается получить в направлениях развития новых методов и алгоритмов решения систем ДАУ с выполнением расчетов на универсальных и на персональных суперкомпьютерах.
Список литературы
1. Библиотека математических программ-решателей на языке Си: SADEL / А.В. Андронов, Д.М. Жук, Д.Ю. Кожевников, В.Б. Маничев. Режим доступа: http://pa10.ru (дата обращения 01.06.2012).
2. Жук Д.М., Маничев В.Б., Андронов А.В. Платформа математического моделирования во временной области разнородных технических систем и объектов
FMS PA10 // Сб. тр. всерос. науч.-тех. конференции «Проблемы разработки перспективных микро и наноэлектронных систем (МЭС)» - 2010. М.: ИППМ РАН, 2010. С. 160-165.
3. Жук Д.М., Маничев В.Б., Ильницкий А.О. Методы и алгоритмы решения дифференциально-алгебраических уравнений для моделирования систем и объектов
во временной области // Информационные технологии. 2010. №7. С. 16-24; № 8. С. 23-26.
4. Жук Д.М., Маничев В.Б. Программа DMAN для решения дифференциально-алгебраических уравнений: свидетельство о государственной регистрации программы для ЭВМ 2009612666; зарегистр. 27.05.2009.
5. Норенков И.П., Трудоношин В.А., Федорук В.Г. Метод формирования математических моделей для адаптируемых программных комплексов анализа радиоэлектронных схем // Радиотехника. 1986. № 9. С.67-72.
6. Маничев В. Б., Глазкова В.Н. Методы интегрирования систем ОДУ для адаптируемых программных комплексов анализа РЭС // Радиотехника. 1988. № 4. С. 88-90.
7. Вержбицкий В.М. Основы численных методов: учебник для вузов. М.: Высшая школа, 2002. 840 с.
8. Решение систем линейных алгебраических уравнений с удвоенной точностью
вычислений на языке Си / В.Б. Маничев, В.Н. Глазкова, Д.Ю. Кожевников, Д.А. Кирьянов, М.К. Сахаров // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. 2011. № 4. С. 25-37.
9. Маничев В.Б. Новые алгоритмы для программ анализа динамики технических
систем // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. 1996. № 1. С. 48-56.
10. Сахаров М. К. Оценка эффективности метода BDF при решении жестких систем обыкновенных дифференциальных уравнений // 14-ая Молодежная международная научно-техническая конференция «Наукоемкие технологии и интеллектуальные системы 2012": труды. М.: МГТУ им. Н.Э. Баумана, 2012. С. 38-41.
11. Numerical Recipes in Fortran: The Art of Scientific Computing / W.H. Press, B.P. Flannery, S.A. Teukolsky, V.T. Vetterling. Cambridge University Press, 1992. 1486 p.
12. Cornea M. Decimal Floating-Point for Intel® Architecture Processors. Режим доступа: http://dx.doi.org/10.1109/ARITH.2009.35 (дата обращения 07.05.2011).
13. Скворцов Л.М. Явный многошаговый метод численного решения жестких дифференциальных уравнений // Журнал вычислительной математики и математической физики. 2007. Т. 47. № 6. С. 959-967.
14. Норенков И.П. Системы автоматизированного проектирования: учеб. пособие для вузов: в 9 кн. Кн. 8. М.: Высшая школа, 1986. 143 с.
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Comparison of modern solvers stiff systems of ordinary differential
equations solvers C library SADEL
# 08, August 2012
DOI: 10.7463/0812.0445558
Zhuk D. M., Manichev V.B., Sakharov M.K.
Russia, Bauman Moscow State Technical University
zhuk@bmstu.ru manichev@bmstu.ru max.sfn90@gmail.com
In this paper the authors present a comparison of various modern solvers of stiff systems of ordinary differentsial equations with C-solvers from SADEL library developed by the authors. A set of test caseswas prepared by the authors; is consists of 6 stiff ODE systems with known analytical or asymptotic solutions. The test results are given in the article. The authors also present the basic theoretical information about the methods used in the testing of the numerical integration of systems of ordinary differential equations and considered various options of problem statement.
Publications with keywords:mathematical modeling, differential algebraic equations (DAE), methods of numerical integration, linear algebraic equations (LAE), ordinary differential equations (ODE)
Publications with words:mathematical modeling, differential algebraic equations (DAE), methods of numerical integration, linear algebraic equations (LAE), ordinary differential equations (ODE)
References
1. Andronov A.V., Zhuk D.M., Kozhevnikov D.Iu., Manichev V.B. Biblioteka atematicheskikh programm-reshatelei na iazyke Si: SADEL [Library of mathematical programs-solvers in C: SADEL]. Available at: http://pa10.ru , accessed 01.06.2012.
2. Zhuk D.M., Manichev V.B., Andronov A.V. Platforma matematicheskogo modelirovaniia vo vremennoi oblasti raznorodnykh tekhnicheskikh sistem i ob"ektov FMS PA10 [Platform of mathematical modeling in the time domain of heterogeneous technical systems and facilities FMS PA10]. Sb. tr. vseros. nauch.-tekh. konferentsii «Problemy razrabotki perspektivnykh mikro i nanoelektronnykh sistem (MES)» - 2010 [Collection of proc. of the all-Russian scientific-technical
conference «Problems of development of prospective micro-and nanoelectronic systems (MES)» -2010]. Moscow, IPPM RAN [Institute of Design Problems in Microelectronics of Russian Academy of Sciences], 2010, pp.160-165.
3. Zhuk D.M., Manichev V.B., Il'nitskii A.O. Metody i algoritmy resheniia differentsial'no-algebraicheskikh uravnenii dlia modelirovaniia sistem i ob"ektov vo vremennoi oblasti [Methods and algorithms of solution of differential-algebraic equations for modeling of systems and objects in the time domain]. Informatsionnye tekhnologii [Information Technology], 2010, no. 7, pp. 16-24; no. 8, pp. 23-26.
4. Zhuk D.M., Manichev V.B. Programma DMAN dlia resheniia differentsial'no-algebraicheskikh uravnenii [DMAN program for solving differential-algebraic equations]. Computer program no. 2009612666, 2009.
5. Norenkov I.P., Trudonoshin V.A., Fedoruk V.G. Metod formirovaniia matematicheskikh modelei dlia adaptiruemykh programmnykh kompleksov analiza radioelektronnykh skhem [Method of forming mathematical models for adaptable software systems of analysis of radioelectronic circuits]. Radiotekhnika, 1986, no. 9, pp. 67-72.
6. Manichev V. B., Glazkova V.N. Metody integrirovaniia sistem ODU dlia
adaptiruemykh programmnykh kompleksov analiza RES [Methods of integration of systems of ordinary differential equations to the adapted program complexes of the analysis of the radio-electronic means]. Radiotekhnika, 1988, no. 4, pp. 88-90.
7. Verzhbitskii V.M. Osnovy chislennykh metodov [Fundamentals of numerical methods]. Moscow, Vysshaia shkola, 2002. 840 p.
8. Manichev V.B., Glazkova V.N., Kozhevnikov D.Iu., Kir'ianov D.A., Sakharov M.K. Reshenie sistem lineinykh algebraicheskikh uravnenii s udvoennoi tochnost'iu vychislenii na iazyke Si [Solving the Linear Algebraic Equation Systems with Double Precision of
Computations in C Language]. VestnikMGTUim. N.E. Baumana. Ser. Priborostroenie [Herald of the Bauman MSTU. Ser. Instrument Engineering], 2011, no. 4, pp. 25-37.
9. Manichev V.B. Novye algoritmy dlia programm analiza dinamiki tekhnicheskikh
System [New algorithms for programs of analysis of the dynamics of technical systems]. Vestnik MGTU im. N.E. Baumana. Ser. Priborostroenie [Herald of the Bauman MSTU. Ser. Instrument Engineering], 1996, no. 1, pp. 48-56.
10. Sakharov M. K. Otsenka effektivnosti metoda BDF pri reshenii zhestkikh sistem obyknovennykh differentsial'nykh uravnenii [Evaluating the effectiveness of BDF method for solving stiff systems of ordinary differential equations]. 14-aiaMolodezhnaia mezhdunarodnaia nauchno-tekhnicheskaia konferentsiia «Naukoemkie tekhnologii i intellektual'nye sistemy 2012": trudy [Proc. of the 14-th international youth scientific and technical conference «Science-Intensive Technologies and Intellectual Systems 2012"]. Moscow, Bauman MSTU Publ., 2012, pp. 38-41.
11. Press W.H., Flannery B.P., Teukolsky S.A., Vetterling V.T. Numerical Recipes in Fortran: The Art of Scientific Computing. Cambridge University Press, 1992. 1486 p.
12. Cornea M. Decimal Floating-Point for Intel® Architecture Processors. Available at: http://dx.doi.org/10.1109/ARITH.2009.35 , accessed 07.05.2011.
13. Skvortsov L.M. Iavnyi mnogoshagovyi metod chislennogo resheniia zhestkikh
differentsial'nykh uravnenii [Explicit multistep method for the numerical solution of stiff differential equations]. Zhurnal vychislitel'noi matematiki i matematicheskoifiziki, 2007, vol. 47, no. 6, pp. 959967. (English version: Computational Mathematics and Mathematical Physics, 2007, vol. 47, no. 6, pp. 915-923. DOI: 10.1134/S0965542507060036 ).
14. Norenkov I.P. Sistemy avtomatizirovannogoproektirovaniia. V9 kn. Kn. 8 [Computer-aided design system. In 9 books. Book 8].Moscow, Vysshaia shkola, 1986. 143 p.
