НАУКА и ОБРАЗОВАНИЕ

УДК.519.622.2

Россия, МГТУ им. Н.Э. Баумана

zhuk@bmstu.ru

manichev@bmstu.ru

max.sfn90@gmail.com

Введение

В математическом ядре любой CAE (ComputerAidedEngineering) системы реализованы различные методы численного интегрирования обыкновенных дифференциальных уравнений (ОДУ) и дифференциально-алгебраических уравнений (ДАУ). Однако, большинство инженеров-проектировщиков не являются специалистами в области численных методов и зачастую не могут оценить адекватность результатов, выдаваемых системой. Поэтому, при компьютерном моделировании и анализе на основе дифференциальных уравнений высокие требования предъявляются к эффективности применяемых численных методов. Эффективность численных методов включает в себя два аспекта:

·       достоверность и точность получаемых решений;

·       затраты машинного времени и оперативной памяти на получение решений.

В данной работе авторы акцентируют внимание на проблеме получения достоверного и точного решения систем ОДУ-ДАУ. Известно, что решение жестких систем ОДУ традиционно вызывает затруднения у большинства численных методов. В наших работах [1-3] было показано, что для гарантии получения корректного и достоверного решения жестких систем ОДУ численный метод их решения должен быть AL-устойчивым, т.е. абсолютно устойчивыми строго в левой (AbsoluteLeft) комплексной полуплоскости. AL-устойчивые методы 2-го и 4-го порядков точности были реализованы в программе DMAN [4].

Главным недостатком реализации вычислительных методов в известных программных комплексах, является возможная выдача ошибочных результатов для жестких систем ОДУ без предупреждения пользователей об их недостоверности. В рамках этой работы жесткой будем считать систему ОДУ со степенью жесткости более 10⁶.

Работа посвящена тестированию различных реализаций методов численного интегрирования, предназначенных для решения жестких систем ОДУ. Предложенный авторами набор тестовых задач был использован для тестирования реализаций методов численного интегрирования в системах MATLAB, Maple, Mathematica, а также в библиотеке NAG и последующего сравнения результатов с решателями библиотеки SADEL (Sets of Algebraic and Differential Equations solvers Library), разработанной авторами [1].

1. Алгоритмы решения систем ДАУ в библиотеке SADEL

В существующих программах-решателях систем ОДУ-ДАУ используются три классические постановки задачи.

1)    Нормальная форма Коши

Где  – вектор дифференцируемых переменных системы размерностью ,  – вектор-функция правых частей,  – независимая переменная (обычно – время)‏. В координатный базис входят только переменные состояния. Такая постановка задачи называется явной. Также заданы начальные условия и отрезок интегрирования , где – заданное значение.

2)    Дифференциально-алгебраическая форма разных индексов

Где  - вектор алгебраических переменных размерностью , - вектор-функция  размерностью . В данном случае в базис входят как дифференцируемые переменные, так и алгебраические. Такой базис называется полным координатным базисом, а постановка задачи – полуявной. Также заданы начальные условия  согласованные начальные условия  и отрезок интегрирования  Начальные условия алгебраических переменных являются согласованными в начальный момент времени, если они удовлетворяют алгебраическим уравнениям.

Отметим, что индекс дифференцирования системы есть наименьшее число аналитических дифференцирований, требующихся для того, чтобы из уравнений (2) получить систему (1). При этом каждое дифференцирование понижает индекс на единицу.

3)    Дифференциально-алгебраическая форма в полном координатном базисе

Где  - вектор-функция размерностью  Заданы начальные условия  согласованные начальные условия  и отрезок интегрирования  Такая постановка задачи называется неявной.

Задача (3) была поставлена в полном координатном базисе Л. Петзольд в 1982 г. Это послужило основой для разработки знаменитой программы DASS (Differential Algebraic Systems Solver) на основе метода BDF (BackwardDifferentialFormula).

Задача решения систем ДАУ в постановке (3), но в расширенном координатном базисе была поставлена в статье [5]. В данном случае в базис помимо дифференцируемых и алгебраических переменных включаются сами производные дифференцируемых переменных. Этот базис используется в решателе систем ДАУ в библиотеке SADEL.

Для решения жестких систем ОДУ были разработаны формулы S-стадийных AL-устойчивых неявных методов интегрирования систем ДАУ вида (3) в расширенном координатном базисе [6]. Формулы имеют вид

где  – число стадий,  – -й шаг интегрирования,  – параметры метода.

На основе этих формул были разработаны и реализованы три метода интегрирования  [4]:

·       М1 – А-устойчивый неявный метод Эйлера первого порядка точности;

·       М2 – АL-устойчивый неявный метод второго порядка точности;

·       М3 – АL-устойчивый неявный метод четвертого порядка точности.

На рисунках 1, 2 проиллюстрирована работа алгоритмов, реализующих эти методы. Методы основаны на совместном решении систем нелинейных алгебраических уравнений  и систем линейных алгебраических уравнений , сформированных на соответствующих стадиях работы алгоритмов [7].

Рис. 1. Схема алгоритмов неявных одно стадийных одношаговых методов Эйлера и трапеций (M1 и M2 - число стадий )

Рис. 2. Схема алгоритма неявного двух стадийного одношагового метода (M3- число стадий )

В библиотеке SADELреализовано три метода решения систем линейных алгебраических уравнений (СЛАУ) [2]. Применение этих методов гарантирует решение с высокой точностью (15 верных значащих цифр) даже плохо обусловленных систем, с числом обусловленности многим больше единицы. Верное решение СЛАУ на каждой итерации алгоритмов интегрирования систем ОДУ-ДАУ способствует верному решение исходной задачи [8].

2. Алгоритмы решения жестких систем ОДУ в современных решателях

Для решения жестких систем ОДУ традиционно используются неявные методы, так как они не имеют ограничения на величину шага интегрирования снизу, при условии, что в программной реализации формул неявного метода  нет арифметических операций деления на шаг [9]. Одним из таких методов является метод BDF.

Фактически, метод BDF – это семейство неявных многошаговых методов численного интегрирования, применяемых для решения жестких систем ОДУ или систем ДАУ [10]. В качестве примера рассмотрим задачу (1). Общая формула для методов BDF может быть записана в виде

где – это шаг интегрирования,  – коэффициент. Величина определяет порядок полинома аппроксимации и число интерполируемых точек до  К примеру,

При  получаем неявный метод Эйлера, при  - метод BDF-2 и т.д. Отметим, что с ростом  область устойчивости метода уменьшается.

Метод BDFреализован практически во всех современных математических средах и библиотеках. В данной работе были протестированы реализации метода в ведущих программных пакетах–WolframMathematica и MATLABSimulink. Помимо этого была протестирована реализация метода в подключаемой библиотеке NAG, написанной на языке Си.

В программном пакете Mapleотсутствует реализация метода BDF. Вместо нее для решения жестких систем ОДУ-ДАУ рекомендовано использование метода Розенброка (Rosenbrock) [11]. Метод Розенброка является L-устойчивым, он обеспечивает хорошую относительную точность  (Tolerance) при небольшой размерности решаемой системы  Рассмотрим задачу (6). Метод Розенброка ищет решение системы в виде

где поправка  находится из решения системы  линейных уравнений. Система имеет вид

где   – матрица Якоби. Коэффициенты  – константы не зависящие от задачи. Если  то получаем метод Рунге-Кутта (Runge-Kutta).

3. Набор тестовых задач ОДУ-ДАУ

Набор тестовых задач состоит из четырех систем ОДУ. Эти системы имеют известное асимптотическое или аналитическое решение. Ниже приведено краткое описание данных тестов.

Тест 1, 1а – это классическая задача Ван дер Поля[12]. Тест представляет собой жесткую систему ОДУ второго порядка, где  – параметр жесткости:

Сравнение программ было проведено для часто встречающихся  на практике значений параметра жесткости (Тест 1) и  (Тест 1а).

Тест 2, 2а – это математическая модель сверхколебательной электронной схемы. Модель представляет собой нежесткую система ОДУ пятого порядка с многопериодным аналитическим решением – тест Маничева [9]:

Здесь

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

Тест 3 представляет собой нелинейную жесткую систему ОДУ второго порядка, имеющую локально неустойчивое решение и составленную по примеру теста Ван дер Поля – тест Скворцова [13]:

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

Тест 4 представляет собой нелинейную жесткую система ОДУ второго порядка для моделирования процессов работы лазера – тест Евстифеева [14]:

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

4. Результаты сравнительного тестирования решателей ОДУ-ДАУ

Тестирование проводилось для параметров решателей, установленных по умолчанию. Все решатели сравнивались при значении относительной точности  В таблице 1 приведены результаты сравнения решателей на предложенном наборе тестовых задач.

Таблица 1. Сводная таблица результатов тестирования решателей ОДУ-ДАУ

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

Необходимо отметить, что реализация метода BDFв WolframMathematicaи реализации в MATLABSimulinkи NAGC-Libraryсправились с разными тестами. Тот факт, что реализация в Mathematicaсправилась с большим количеством тестов, свидетельствует о модификации метода BDF разработчиками Wolfram.

При дополнительных настройках параметров решателей удалось существенно улучшить результаты тестирования. Например, при искусственном увеличение точности решатель библиотеки NAGне справился лишь с тестом 2а, а техника комбинирования явного метода Рунге-Кутта и неявного метода Эйлера при повышенной точности в пакете WolframMathematicaпозволила успешно решить весь набор тестовых задач. Однако такие настройки требуют от инженеров-проектировщиков широких знаний в области методов численного интегрирования систем ОДУ и хорошего владения используемыми математическими пакетами. Так как серьезной проблемой является почти полное отсутствие диагностических сообщений о возможных ошибках, то инженер должен также уметь оценить адекватность полученного решения. Решатели библиотеки SADELпозволяют решать любые системы ОДУ-ДАУ, выдавая гарантировано верный результат в широком диапазоне значений параметров жесткости, что существенно облегчает жизнь простому инженеру.

На рисунках 3 – 6 представлены некоторые результаты сравнения решателей ОДУ-ДАУ.

Рис. 3. Решение теста 2а полученное с помощью библиотеки NAGпри

Рис. 4. Правильное решение теста 2а

Рис. 5. Решение теста 3 полученное с помощью библиотеки NAG при

Рис. 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 с.