Библиотека решателей

Перейдем к рассмотрению специальных алгоритмов численного анализа непрерывных режимов, сконструированных с учетом особенностей

обозначенных классов систем (1.4) и (1.5). При анализе жестких задач в современных инструментальных средах [54, 65, 74] обычно рекомендуется применение неявных формул (неявный метод Рунге-Кутта, формулы Розенброка, алгоритмы RADAU5 и DDASSL).

Для решения задач (1.4) будем использовать одношаговые схемы, которые имеют вид

где a_i,- числовые коэффициенты, k_i- стадии

метода, m- число стадий. Численные формулы (2.5) хорошо приспособлены для практических расчетов, поскольку они не требуют вычисления дополнительных начальных значений и позволяют легко изменять шаг интегрирования.

Многие современные алгоритмы на основе явных методов с контролем точности не приспособлены для решения жестких задач, поскольку такой подход приводит к потере эффективности и надежности из-за «раскачивания» шага интегрирования на участке установления вследствие потери

74 устойчивости.

В работах [57, 70] рассматривается абсолютная устойчивость метода. Понятие абсолютной устойчивости вводится на линейном скалярном уравнении с комплексным.

Точное решение задачи (2.6) в точке t_n+1 записывается в виде

при этом выполняется неравенствоВ

результате применения метода (2.5) для решения (2.6) получается

- функция роста или функция устойчивости. Функция устойчивости m-стадийного метода (2.5) есть многочлен степени m [57, 76]:

Очевидно, что при численном решении (2.6) также должно быть истинным неравенствоМетод, удовлетворяющий этому условию для

данного z, называется абсолютно устойчивым. Область Rкомплексной плоскости Zназывается областью устойчивости метода, если он устойчив при всех z ∈ R. В дальнейшем будем использовать понятие интервала устойчивости - пересечение области устойчивости с вещественной осью [57] комплексной плоскости Z.

Поскольку при z ®-∞ имеет местообласти устойчивости

явных методов типа Рунге-Кутты ограничены.

постоянная D- длина интервала устойчивости, а l_max- максимальное собственное число матрицы Якоби задачи (1.4).

Как правило, процесс исследования устойчивости является достаточно трудоемким. Для автоматизации решения данной задачи, на основе подхода, описанного в [15, 54, 57], автором разработаны алгоритмы и программные средства [17, 73] анализа устойчивости явных численных схем, которые обеспечивают решение следующих задач: 1) расчет коэффициентов полинома устойчивости по заданным параметрам численной схемы; 2) построение области устойчивости заданной конфигурации, поиск соответствующей функции роста и расчет параметров численного метода с заданной областью устойчивости; 3) графическая интерпретация функции устойчивости

(построение области устойчивости); 4) расчет длины интервала устойчивости; 5) создание отчетов об исследовании в формате Microsoft Word или Excel.

Комплекс программ исследования устойчивости интегрирован в среду моделирования динамических и гибридных систем ИСМА [73] и зарегистрирован в Роспатент [17]. Правильность работы программ доказана конструктивно на всем множестве тестовых примеров из [57]. Вместе с малостадийными методами разного порядка точности успешно исследованы многостадийные схемы Фельберга, в том числе с расширенной областью устойчивости [53, 56]. Использование созданного инструментария позволяет существенно сократить временные затраты исследователя, занимающегося созданием эффективных алгоритмов численного анализа на основе явных схем.

Приведем краткую характеристику (таблица 2.2) алгоритмов интегрирования, включенных в библиотеку решателей ИСМА [19]. Наряду с классическими широко известными численными схемами используются и оригинальные алгоритмы. Это явные схемы переменного порядка и шага с контролем точности и устойчивости [53, 54, 55, 56, 57, 74, 88, 96], неявные методы, а также полуявные (m-k)-схемы второго порядка точности с замораживанием матрицы Якоби [52, 58, 81, 93].

точности, сочетающий многостадийные явные схемы и неявный метод

RADAU5.

Таблица 2.2 - Решатели среды ИСМА

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

виде API. Используя API, любой пользователь с начальными знаниями объектно-ориентированного программирования имеет возможность разрабатывать и встраивать в систему новые численные методы без перекомпиляции всей системы в целом.

Алгоритм DISPSl

Рассмотрим подробнее особенности алгоритма DISPS1, авторская реализация которого включена в набор решателей ИСМА [16, 25]. В работах [28, 54, 55] показано, что повысить эффективность алгоритма интегрирования можно за счет использования численных схем (2.5) с различными областями устойчивости. Так, на участке установления, когда критическим фактором является устойчивость метода, необходимо использовать схему с возможно большим интервалом устойчивости. С другой стороны, на участке возмущения, где шаг интегрирования определяется точностью, целесообразно перейти к более точной схеме. Для увеличения интервала устойчивости следует использовать многостадийные численные схемы. Однако с ростом числа стадий пропорционально увеличиваются вычислительные затраты, связанные с вычислением правой части решаемой задачи. Поэтому применение на участках установления столь же точных методов, как на участке возмущения, но с большим числом стадий, является неэффективным. Но, как показано в [57], на участке установления можно использовать методы меньшего порядка точности, при этом глобальная ошибка решения не ухудшается. Объясняется это тем, что, как правило, большая часть интервала решения связана с участками установления, где шаг выбирается из условия устойчивости, в результате чего точность оказывается лучше требуемой. Таким образом, интегрирование на участке установления целесообразнее производить методами низких порядков точности.

В [57] предложен алгоритм переменного порядка и шага с контролем точности и устойчивости DISPS, высокая эффективность которого достигается применением наиболее выгодной численной схемы в зависимости от поведения

решения в процессе интегрирования. Для этого сформулированы требования к параметрам используемых методов и построены условия для контроля точности и устойчивости.

Пусть заданы наборы параметров численных схем первого, второго и третьего порядков точности с различным числом стадий. В отличие от алгоритма, приведенного в [57], допускается, что число стадий в двух соседних методах одного порядка может отличаться на любую величину.

Для хранения используемых схем введем в рассмотрение матрицу методов размерностью MaxK ? MaxM, где MaxK- максимальный порядок точности, MaxM- максимальное число стадий по всем используемым методам. Поскольку минимальное число стадий ограничено тремя, первые два столбца матрицы не используются, и при программной реализации их можно не задавать. Элемент матрицы методов на пересечении к -й строки и m-го столбца содержит коэффициенты и интервалы устойчивости для m-стадийного метода к -го порядка. На рисунке 2.16 приведена таблица, соответствующая матрице методов. В ее ячейках записаны интервалы устойчивости используемых методов.

Рисунок 2.16 - Матрица методов

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

Поскольку алгоритм DISPS1 может включать различное число методов любого порядка, целесообразно предоставить пользователю возможность самому определять набор используемых методов. Для этого необходимо, используя модуль исследования устойчивости, заранее сформировать множество схем с желаемыми областями устойчивости. Затем из готовых схем в конфигураторе DISPS1 набирается конкретная реализация алгоритма. Число таких реализаций может быть любым. При этом добавление нового метода или набора методов осуществляется без перекомпиляции. При загрузке библиотеки алгоритмов, содержащей реализацию DISPS1, сформированные ранее наборы становятся доступными для использования. Таким образом, пользователь имеет в распоряжении библиотеку алгоритмов интегрирования, содержащую набор готовых методов, и может самостоятельно включать в нее собственные численные схемы.

Результаты тестирования алгоритма с адаптивной областью устойчивости, приведенные в работах [16, 22, 57, 82], служат конструктивным доказательством высокой эффективности решателя. Следует отметить хорошие результаты решения всего множества многомерных задач из [94]. Поэтому метод DISPS1, как и другие решатели переменного порядка среды ИСМА, может быть успешно использован для анализа переходных процессов в ЭЭС, описываемых системами ОДУ вида (1.4).