авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 |

«Министерство образования и науки Российской Федерации САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ Ю.Б. Сениченков ...»

-- [ Страница 3 ] --

Иначе 3. [t i*1, ti* ] : Длительное поведение. Решение уравнений ds k = f k (t, s k, xk ) dt 0 = g k (t, s k, xk, y k ) x1 = y 2 ;

x2 = y k = 1, от согласованных начальных условий, до тех пор, пока не станет истинным хотя бы один предикат pred k ( s k (t ), xk (t ), y k (t )).

Аварийный выход: не смогли найти решить уравнения.

Как только хотя бы один предикат стал истинным, выполняем:

завершается старый интервал гибридного времени i t i* := t : pred k ( s (t ), x(t ), y (t )) = true pre( s k (t i* )) := s k (t );

pre( xk (t i* )) := xk (t );

pre( y k (t i* )) := y k (t );

начинается новый интервал гибридного времени i + Переходим к 2.

4. Post _ gapi : Завершающие мгновенные действия.

Если a) pred1 = true и pred 2 = false то post ( s1 (t i*1 )) := ( Init1 ( pre( s1 (t i*1 )), pre( x1 (t i*1 )), pre( y1 (t i*1 )) post ( s 2 (t i*1 )) := pre( s 2 (t i*1 )) переходим к d) здесь post ( y k (t i*1 )) - новое начальное приближение для поиска новых начальных условий.

Если b) pred1 = false и pred 2 = true то post ( s1 (t i*1 )) := pre( s1 (t i*1 )) post ( s 2 (t i*1 )) := Init 2 ( pre( s 2 (t i*1 )), pre( x 2 (t i*1 )), pre( y 2 (t i*1 )) переходим к d) d) post ( xk (t i*1 )) := pre( xk (t i*1 ));

post ( y k (t i*1 )) := pre( y k (t i*1 )) k = 1, Переходим к 2.

Конец алгоритма_D.

Если бы формулы (5) вычислялись бы одновременно, то есть сначала бы вычислялись новые значения начальных условий на всех одновременно сработавших переходах, и только потом новые значения связей, то мы бы вновь получили периодическое движение по обеим компонентам. Именно так реализована обработка событий, приводящих к смене поведения в последней версии Model Vision Studium.

Непрерывная стратегия синхронизации может быть представлена следующим образом.

Алгоритм_С работы функциональной схемы из двух примитивных автоматов с непрерывной синхронизацией одновременных событий Начало Алгоритма_C 1. Инициализация.

Начинается первый промежуток t0 = * post ( sk (t 0 )) := pre( s k (t 0 )) := sk 0 ;

* * post ( xk (t 0 )) := pre( xk (t 0 )) := xk 0 ;

* * post ( y k (t 0 )) := pre( y k (t 0 )) := y k * * k = 1, здесь - y k 0 или согласованное начальное условие, либо только начальное приближение к нему.

1. Pr e _ gapi : Вычисление новых согласованных начальных условий и проверка предиката.

i := i + 1;

Решаем уравнения 0 = g k (t i*1, post ( s k (t i*1 )), post ( xk (t i*1 )), post ( y k (t i*1 ))) k = 1,2;

дополненные уравнениями связи post ( x1 (t i*1 )) = post ( y 2 (t i*1 ));

post ( x 2 (t i*1 )) = post ( y1 (t i*1 ));

относительно post ( y k (t i*1 )) с начальным приближением pre( y k (t i*1 )).

Аварийный выход: не смогли найти согласованные начальные условия.

Вычисляем предикаты (одновременно) на левом конце нового промежутка pred k ( post ( s (t i*1 )), post ( x (t i*1 )), post ( y (t i*1 )) Если хотя бы один предикат истинен pred k ( post ( s(t i*1 )), post ( x (t i*1 )), post ( y (t i*1 )) = true, то завершается текущий интервал гибридного времени i t i* := t i* pre( s k (t i* )) := post ( s k (t i* ));

pre( xk (t i* )) := post ( xk (t i* ));

pre( y k (t i* )) := post ( y k (t i* ));

переходим к 4.

Иначе 3. [t i*1, t i* ] : Длительное поведение. Решение уравнений ds k = f k (t, s k, xk ) dt 0 = g k (t, s k, xk, y k ) x1 = y 2 ;

x2 = y k = 1, от согласованных начальных условий, до тех пор, пока не станет истинным предикат хотя бы один предикат pred k ( s k (t ), xk (t ), y k (t )).

Аварийный выход: не смогли найти решить уравнения.

Как только хотя бы один предикат стал истинным, вычисляем все предикаты (одновременно) на правом конце промежутка и выполняем завершается старый интервал гибридного времени i t i* := t : pred k ( s (t ), x(t ), y (t )) = true pre( s k (t i* )) := s k (t );

pre( xk (t i* )) := xk (t );

pre( y k (t i* )) := y k (t );

начинается новый интервал гибридного времени i + Переходим к 2.

4. gapi : Завершающее мгновенное поведение.

a) pred1 = true и pred 2 = false post ( s1 (t i*1 )) := ( Init1 ( pre( s1 (t i*1 )), pre( x1 (t i*1 )), pre( y1 (t i*1 )) post ( s 2 (t i*1 )) := pre( s 2 (t i*1 )) переходим к d) здесь post ( y k (t i*1 )) - новое начальное приближение для поиска новых начальных условий.

b) pred1 = false и pred 2 = true post ( s1 (ti*1 )) := pre( s1 (ti*1 )) post ( s 2 (t i*1 )) := Init 2 ( pre( s 2 (t i*1 )), pre( x 2 (t i*1 )), pre( y 2 (t i*1 )) переходим к d) c) pred1 = true и pred 2 = true post ( s1 (t i*1 )) := ( Init1 ( pre( s1 (t i*1 )), pre( x1 (t i*1 )), pre( y1 (t i*1 )) post ( s 2 (t i*1 )) := ( Init 2 ( pre( s1 (t i*1 )), pre( x 2 (t i*1 )), pre( y 2 (t i*1 )) переходим к d) d) post ( xk (t i*1 )) := pre( xk (t i*1 ));

post ( y k (t i*1 )) := pre( y k (t i*1 )) k = 1, Переходим к 2.

Конец алгоритма_C.

Очевидно, что выбор частного случая блок-схемы в виде кольца, не влияет на работу обоих алгоритмов для произвольных схем. Все отличия связаны только с выбором первого обрабатываемого события для первого, и проверкой всех переходов, исходящих из текущего узла – во втором.

2.11 Открытый гибридный автомат с контактами.

Компонентные модели можно строить, используя не только «входы выходы», но и контакты. В изолированном примитивном гибридном автомате H = {w, f ( w), g ( w),W0, pred (t, w(t )), Inv(t, w(t )), Init ( w(t ))}, w = {s1, s 2 } с системой уравнений ds f ( s1, s 2, t ) = 0;

dt g ( s1, s 2, t ) = где s1 «дифференциальные» и s2 «алгебраические» составляющие вектора переменных состояний w, представим множество переменных w в виде двух непересекающихся множеств w = Cont S и Cont S =, и будем называть элементы z множества Cont контактами, а элементы s множества S - переменными состояния. Точно так же как в автомате с входами и выходами, переменные состояния s автомата с контактами инкапсулированы в нем и не могут изменяться извне. Связывание автоматов в блок-схемы происходит с помощью контактов z. В отличие от входов и выходов, место которых строго фиксировано в структуре уравнений - а именно входы и выходы могут входить только в правые части дифференциальных уравнений и формул, контактом может быть названа любая переменная состояния автомата.

Hс1O O Два открытых примитивных автомата с котактами и Hс Hс1O Hс2, если хотя бы O называются объединенными в блок-схему контакт одного их них равен контакту другого (связаны между собой).

Соответствующие пары переменных называются связными. Не связные переменные называются свободными.

Композицией Hс10 Hс2 двух примитивных автоматов Hс1O и Hс1O с контактами z1 и z 2 называется автомат с новой системой уравнений dz1d f11 ( w1, t ) = 0;

dt ds1d f12 ( w1, t ) = 0;

dt g1 ( w1, t ) = dz 2 d f 21 ( w2, t ) = 0;

dt ds2 d f 22 ( w2, t ) = 0;

dt g 2 ( w2, t ) = дополненной уравнениями контактов:

z1c = z 2 c ;

z1c z1 ;

z 2 c z 2, с новым предикатом, представляющим объединение предикатов pred 1 (t, w1 (t )) or pred 2 (t, w1 (t )) и новой функцией инициализации, где через s1d, s2 d, z1d, z 2 d обозначены дифференциальные компоненты векторов состояния. Выделение в отдельную группу дифференциальных уравнений относительно неизвестных контактов обусловлено только тем, что при этом возникают новые задачи при построении совокупной системы.

Прежде всего, при построении итоговой системы мы получаем готовую к реализации модель, ее уравнения можно передавать непосредственно решателям. Блоки с контактами имеют другое назначение – это компоненты схем, и итоговая система обычно должна быть предварительно преобразована. Исходная система для компоненты обычно оказывается переопределенной. В отдельных блоках может полностью отсутствовать информация о том, какие переменные считать неизвестными. Например, рассмотрим блок с контактами z1 и z 2, и dz уравнением z1 = 0. В нем отсутствует информация о том, следует ли dt dz считать написанное уравнение дифференциальным - = z1 относительно dt dz переменной z 2 и правой частью z1 или формулой - z1 =, в которой dt переменная z1 является производной от некоторой функции z 2.

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

2.12 Частичная синхронизация гибридных автоматов.

В предыдущих параграфах рассматривались композиция автоматов H 10 и H 2, подразумевавших полную их синхронизацию в гибридном времени. Однако во многих случаях полная синхронизация является слишком жестким требованием – реальные объекты и их модели могут обмениваться информацией только в определенные моменты своего гибридного времени. И только в эти моменты, моменты частичной синхронизации и могут влиять друг на друга гибридные автоматы. В тоже время между автоматами может не происходить обмена информации вовсе. Поэтому целесообразно ввести еще одно понятие – совместно работающую пару частично синхронизированных автоматов H 1O H 2, O подразумевая под этим, что каждый автомат работает в своем гибридном времени 1 и 2, и нет возможности установить одновременность любых двух событий в этих автоматах. Это крайний случай интересен тем, что решаемые системы уравнений не связаны между собой, и их можно решать параллельно различными численными методами в своем модельном времени.

Для того чтобы ввести понятие частичной синхронизации, рассмотрим 1 = {( pre_ gapk,[tk, tk +1 ], gapk )} последовательности и 2 = {( pre _ gapj,[t j, t j +1 ], gapj )} двух автоматов. Назовем синхронизирующей последовательностью последовательность пар gapk s = {(gapk 1, gapj 2 ) m,}, gap j где и являются упорядоченными по возрастанию подпоследовательностями в своем гибридном времени.

Назовем два примитивных изолированных или открытых гибридных автомата H 10 и H 2 парой частично синхронизированных автоматов, понимая под этим автомат, определенный на дискретном времени T = {0,1,2...} с вектором переменных ws = [ w1, w2 ] и неявно заданным F : wsm+1 = F ( wsm ), определяемым уравнениями каждого отображением отдельного гибридного автомата на своих подпоследовательностях gap k 1 и gap j 2.

2.13 О практической ценности примитивных автоматов.

Примитивные гибридные автоматы могли бы представлять лишь теоретическую ценность, если бы в программном комплексе требовалось бы сводить автоматы со многими состояниями и образованные из них блок-схемы к примитивным автоматам. На счастье этого делать не надо, так как сама природа гибридных автоматов подсказывает иное, очевидным образом более экономичное решение.

Представленные выше автоматы с несколькими состояниями представляют последовательные процессы – в каждый момент времени текущим является только один узел. Следовательно, текущий узел и принадлежащие ему выходные переходы полностью определяют поведение гибридного автомата в каждый момент. Таким образом, в каждый момент автомат с несколькими состояниями эквивалентен примитивному автомату. Если учесть, что в детерминированном автомате с нескольким узлами может одновременно сработать только один переход, то мы получаем «динамическую» модель, для работы которой необходима только композиция примитивных гибридных автоматов. При этом нет необходимости анализировать структуру графа и сводить его к трудно проверяемой теоретической схеме.

Литература 1. Колесов Ю.Б., Расторгуев И.В., Сениченков Ю.Б. Современное программное обеспечение для проведения вычислительного эксперимента./ В сб. «Труды Технического университета», № 468, 1997, стр. 95-98.

2. Yu. Kolesov, Yu. Senichenkov. Model Vision 3.0 for Windows 95/NT. The graphical environment for complex dynamic system design. ICI&C’97.

International conference on Informatics and Control. St. Petersburg, 1997. pp.

764-768.

3. Yu. Kolesov, Yu. Senichenkov. Visual Specification language intended for event-driven hierarchical dynamical systems with variable structure ICI&C’97.

International conference on Informatics and Control. St. Petersburg, 1997. pp.

704-711.

4. Yu. Kolesov, Yu. Senichenkov. Model Vision 3.0 for WINDOWS 95/NT - graphical environment for complex dynamic systems modeling. / COLOS’98.

Maribor 1998, pp.34-41.

5. Колесов Ю.Б., Сениченков Ю.Б. Семейство Model Vision.

«Информационные технологии в моделировании и управлении». II Межд.

Конференция, С.-Петербург, 2000, стр.188- 6. Колесов Ю.Б., Сениченков Ю.Б. Современные подходы к компьютерному моделированию динамических систем. Научно технические ведомости СПбГТУ, 3' 2002. с.93-102.

7. Yu. Kolesov, Yu. Senichenkov. A composition of open hybrid automata.

EUROCON 2003. IEEE Catalog Number: 03EX665, ISBN: 0-7803-7763-X, Library of Congress: 2002117046.

8. Колесов Ю.Б., Сениченков Ю.Б. Компьютерное моделирование в научных исследованиях и в образовании, Exponenta_Pro. Математика в приложениях. 1' 2003. с.4- 9. Колесов Ю.Б., Сениченков Ю.Б. Синхронизация событий при использовании гибридных автоматов для численного моделирования сложных динамических систем. Научно-Технические ведомости, 1’ 2004, стр. 202- Глава 3. Формирование итоговой системы уравнений и продвижение модельного времени.

Состав библиотеки численных методов существенно зависит от типа и формы уравнений, передаваемых решателю, требований, предъявляемых к численному решению и набора инструментов, используемых при исследовании свойств модели.

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

В тоже время от формы и типа уравнений во многом зависит эффективность решения. Анализом формы итоговой системы в MVS занимается блок структурного анализа уравнений.

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

Набор инструментов для исследования свойств модели определяется входным языком MVL. Эти инструменты либо входят в базовое программное обеспечение, либо оформляются пользователем в виде специальных библиотек (DLL).

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

Первыми типами блоков, использовавшимися в семействе Model Vision, были блоки «вход-выход-состояние». Принятой формой описания поведения такого блока является система уравнений ds = f ( x, s, t ), dt y = g ( s, x, t ) где s - вектор переменных состояния, x – входной, а y – выходной векторы.

В ранних версиях пакета предполагалось, что пользователь пишет уравнения именно в такой форме, приемлемой для любого численного метода, из числа существующих в библиотеке MVS. Однако соединение таких блоков в функциональные схемы может приводить к алгебро дифференциальным уравнениям ds = f ( x, s, t ), dt 0 = g ( y, s, x, t ) что заставило считать последнюю форму также допустимой. При использовании программ, реализующих явные методы, она требует всего лишь дополнительно программных реализаций методов решения нелинейных уравнений и может быть непосредственно передана решателю при использовании программных реализаций неявных методов.

Появившиеся позже блоки с контактами привели к свободной форме ds, s, t ) = 0, F( dt в которой уравнения первого порядка не разрешены относительно первых производных.

Очевидно, что системы уравнений более высокого порядка d ns ds F ( n,...,, s, t ) = dt dt приводятся к свободной форме автоматически. Ориентация на свободную форму оправдана еще и тем, что это не только удобная форма представления исходной информации, но и форма, естественным образом возникающая при «физическом моделировании» [1].

Свободная форма редко используется в программных реализациях численных методов [2]. Более распространенной является форма ds = f ( s, t ), A( s ) dt с возможно вырожденной матрицей A(s ).

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

Структурному анализу можно подвергать как системы отдельного блока, так и всю, совокупную систему уравнений. На уровне отдельного блока наиболее часто возникает следующая ситуация. Разработка блока начинается с создания класса «Система уравнений», в котором переменные еще не делятся на входы, выходы, контакты и внутренние переменные.

Система уравнений может оказаться недоопределенной или переопределенной, и можно либо потребовать от пользователя уточнения списка искомых переменных, после чего проверить ее разрешимость, либо можно попытаться построить все возможные разрешимые системы, и показать их пользователю. На уровне отдельного блока последняя постановка вполне реалистична, так как обычно число переменных и уравнений отдельного блока не очень велико. Такой подход принят в пакете ABACUSS [3].

Сформированная система передается решателю, задача которого найти численное решение в конкретной временной точке. Выбор этой конкретной точки определяется блоком продвижения модельного времени и зависит от режима работы пакета.

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

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

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

Рассмотрим работу отдельных блоков подробнее.

3.1 Структурный анализ уравнений.

Рассмотренная во второй главе форма представления уравнений примитивных автоматов чрезвычайно сложна для восприятия человеком.

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

Альтернативная форма представления в виде автомата со многими длительными состояниями проще для восприятия. Это объясняется тем, что каждому узлу приписана обычно легко узнаваемая система уравнений небольшой размерности, а структура графа отражает внутреннюю логику поведения изучаемого объекта. Одного этого достаточно, чтобы отказаться от использования примитивных автоматов в качестве основной конструкции графического языка. Однако есть еще более важное обстоятельство – как уже отмечалось, при использовании гибридных автоматов со многими состояниями в каждый момент смены состояния заново формируется, а затем решается система небольшой размерности без дополнительных вспомогательных уравнений, что приносит значительный вычислительный эффект.

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

Процедура формирования совокупной системы и ответ на вопрос о ее разрешимости зависят от типа используемых блоков. В MVS возникают следующие типы блок-схем:

• несвязанные между собой блоки с входами-выходами или контактами, • блок-схемы из связанных блоков с входами-выходами, • блок-схемы из связанных блоков с контактами, • блок схемы из блоков произвольного типа.

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

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

В работе [4] приводятся наиболее распространенные формы записи и канонические формы для систем алгебро-дифференциальных уравнений.

Первый и наиболее простой тип систем - это системы линейных уравнений ds + B * s = f (t ) с постоянными коэффициентами. Будем считать, что A* dt матрицы входящие в эти уравнения являются вещественными квадратными матрицами. Это форма называется неявной или формой с неразрешенными относительно производных линейными уравнениями.

Более сложным являются линейные уравнения с переменными ds + B (t ) * s = f (t ). В обоих случаях вводят еще и коэффициентами A(t ) * dt полуявную форму ds + B11 (t ) s1 + B12 (t ) s 2 = f dt B21 (t ) s1 + B22 (t ) s 2 = f выделяя в векторе неизвестных s дифференциальные s1 и алгебраические s2 переменные. Нелинейные системы обычно представляют в свободной ds ds, s, t ) = 0 и ее модификации + g ( s ) = f (t ),когда для форме F( A* dt dt разрешения системы относительно производных требуется решение системы линейных алгебраических уравнений.

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

Рассмотрим две системы, состоящие из двух уравнений:

dx x y = dt 0 = x + 2 y + f (t ) и dx dy + f1 (t ) = dt dt 0 = x + y 2 f 2 (t ) Первая из них может быть непосредственно передана решателю, а вторая требует дифференцирования второго уравнения dx dy d + 2 y + f 2 (t ) = 0.

dt dt dt dx После чего, подставляя значение из исходного дифференциального dt dy dy d + f1 (t ) + 2 y + f 2 (t ) = уравнения в новое, находим или dt dt dt dy d (2 y 1) = f1 (t ) f 2 (t ). В обеих формах речь идет об одном dt dt «независимом» дифференциальном уравнении и одной функциональной зависимости, но во втором случае требуется анализ и дополнительное предварительное преобразование.

Еще более ярким является пример системы [5] dx dy + x = f (t ) dt dt dx dy + +y= dt dt с решением d x(t ) = f (t ) f (t ) dt.

d y (t ) = f (t ) dt В данном случае исходная система из двух дифференциальных уравнений на самом деле оказывается уравнением и функциональной зависимостью.

Уравнения dx + y = f1 (t ) dt x = f 2 (t ) вообще не содержат ни одного дифференциального уравнения.

При «физическом» подходе к моделированию, системы, аналогичные рассмотренным выше, появляются практически всегда [6].

Рассмотрим приводимый авторами Modelica пример. Пусть мы имеем два готовых устройства — электродвигатель постоянного тока M и пропеллер P — и хотим построить из них систему-вентилятор, соединив их жестким валом (рис. 1) Рис. 1. Соединение мотора и пропеллера в единую систему-вентилятор В данном примере мы имеем два очевидных блока — двигатель M и пропеллер P, поведение которых упрощенно задается следующими уравнениями:

d k1 I + T dt = J d T k p a = (M) (P), dI = V k 2 Ra I dt Jp dt La где, T - угловая скорость и момент на валу двигателя и пропеллера, V — напряжение на двигателе, I — ток якоря, J a, Ra, La — момент инерции, омическое и индуктивное сопротивление якоря соответственно, J p момент инерции пропеллера, k1, k 2, k p - специальные коэффициенты.

Очевидно, что переменные M.V, M., M.T являются внешними для блока M, а переменные P., P.T — для блока P, в том смысле, что мы имеем модели «мотора, вращающего нагрузку» и «нагрузки, вращаемой мотором». И если задать коэффициенты дифференциальных уравнений мотора и нагрузки, то мы получим классические динамические системы независимых физических объектов.

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

M. = P.

M.T + P.T = Полученная итоговая система является алгебро-дифференциальной и требует дифференцирования одного из уравнений.

Одним из возможных способов решения систем алгебро дифференциальных уравнений рассмотренного типа является приведение исходной системы к новой системе дифференциальных уравнений дифференцированием алгебраических уравнений. Число необходимых дифференцирований, после которых полученная новая система будет системой дифференциальных уравнений относительно всех входящих в нее неизвестных, называется индексом системы. Например, для системы dx1 dx + = sin(T ) dt dt x1 + 2 x 2 = x1 (0) = 1;

x2 (0) = x1 (0);

индекса 1, новой системой будет dx1 dx + = sin(T ) dt dt dx1 dx +2 2 = dt 1 dt x1 (0) = 1;

x2 (0) = x1 (0);

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

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

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

2. Найти ее индекс и построить систему, передаваемую решателю.

3. Найти согласованные начальные условия для построенной системы.

В зависимости от типа системы:

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

Для линейных систем с постоянными коэффициентами существует каноническая форма, позволяющая одновременно и определить индекс системы и построить ее решение. Любая система это типа может быть приведена к форме dy + PBQy = f PAQ dt E 0 C PAQ = ;

PBQ = 0 N 0 E где P, Q - квадратные невырожденные матрицы, E - единичная, N нильпотентная матрица, индекса k. Это позволяет представить исходную систему в виде dy + Cy1 = f dt Ny 2 + y 2 = f где первая система – это система обыкновенных дифференциальных k y 2 = (1) i N i f 2(i ), где f 2(i ) уравнений, а вторая – набор формул i = производные вектора правой части.

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

Рассмотрим систему алгебро-дифференциальных уравнениях вида dx + Bx = f, x, f n, A, B n,n A dt с постоянными коэффициентами и будем рассматривать только квадратные матрицы, так как случай прямоугольных матриц не вносит никаких особенностей.

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

Рассмотрим случай, когда ранг матрицы отличен от нуля и меньше числа dxi vi = ;

C = [ A, B];

z = [v, x]T и будем уравнений. Введем обозначения dt приводить к треугольному виду матрицу C системы Cz = f.

Так как ранг матрицы C меньше числа уравнений, то она примет следующий вид показанный на Рис.2. Обозначим матрицы A и B и вектор f после преобразования через A1, B1, f1 соответственно. Если матрица B2, состоящая из последних m строк матрицы B1, имеет ранг m, то систему B2 x = f m можно разрешить относительно m компонентов вектора x и подставить найденные выражения в полученную систему dx + B11 x = f n1m, где A11, B1 - матрицы, содержащие первые n-m строк A11 dt матриц A1, B1, а вектор f n1m - первые n-m компонентов вектора f1.

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

Однако, если матрица B2 не является матрицей полного ранга, то здесь возможны два случая:

• При приведении матрицы B2 к трапециевидной появляются нулевые строки и соответствующие им нулевые компоненты вектора правой части. В этом случае система не доопределена.

• При приведении матрицы B2 к трапециевидной появляются нулевые строки и соответствующие им не нулевые компоненты вектора правой части. В этом случае система не совместна.

A B n n m 0 m Рис 2. Вид матрицы С после первого преобразования у треугольному виду.

Таким образом, за конечное число шагов либо будет построена система, имеющая единственное решение, дополненная системой формул, либо показано, что систему нельзя передавать решателю.

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

К сожалению, в общем случае можно проанализировать только структуру уравнений, ответив на вопрос, не является ли написанная система системой структурного ранга больше единица. При этом используется только информация о вхождении переменных в то или иное уравнение, но не собственно коэффициенты, что приводит к необходимости говорить о структурном ранге. Примером подобного анализа является предварительный анализ структурной вырожденности больших систем линейных алгебраических уравнений [8]. Напомним, что система линейных алгебраических уравнений с разреженной квадратной Ax = b называется структурно вырожденной, если ее матрицей A:

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

Структура системы определяется матрицей инцидентности, строки которой соответствуют уравнениям системы, столбцы – неизвестным.

Элементы матрицы имеют нулевые значения в соответствующем столбце, если данная переменная не входит в рассматриваемое уравнение, и единицы – в противном случае.

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

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

ds, s, t ) = 0 могут быть переданы в Уравнения в свободной форме F ( dt MVS только программе DASSL. Уравнения в полуявной форме (дифференциальные уравнения с ограничениями или система индекса 1 в форме Хессенберга) ds = f (t, s, y ) dt 0 = g (t, s, y ) g при невырожденной матрице, могут быть переданы любому z численному методу, из числа имеющихся в библиотеке.

Часто рассматривают также практически важный случай систем индекса в форме Хессенберга ds = f (t, s, y ), dt 0 = g (t, s ) однако эта система уже требует дифференцирования.

Полуявная форма может также быть представлена в виде (метод фиктивных производных) [8] dw =y dt ds dw = f (t, s, ).

dt dt dw 0 = g (t, s, ) dt Для того чтобы интегрировать уравнения в свободной форме любым численным методом, их преобразуют к полуявной форме [8] ds =z, dt 0 = f (t, s, z ) при этом индекс увеличивается на единицу по отношению к исходной системе.

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

Constantios C. Pantelides. The consistent initialization of differential-algebraic systems. SIAM J. Sci. Stat. Comput. Vol.9, № 2, March 1988.

В этой работе ставится вопрос о том, какое минимальное число уравнений исходной системы надо продифференцировать, чтобы найти независимые дифференциальные уравнения и функциональные зависимости, и предлагается следующая форма исходных уравнений, которая для системы dx1 dx + = sin(T );

dt dt x1 + 2 x2 = 0;

x1 (0) = 1;

x2 (0) = x1 (0);

будет выглядеть так dx2 dx z2 = ;

z1 = 1 ;

z1 + z 2 = sin(T );

dt dt x1 + 2 x2 = 0;

x1 (0) = 1;

x2 (0) = x1 (0) Эта система после анализа и переформирования превратиться в систему, dx2 dx z2 = = z1 ;

z1 = sin(T ) z 2 ;

;

dt dt, x1 + 2 x2 = 0;

x1 (0) = 1;

x2 (0) = x1 (0);

dx требующую либо явного выражения только для производной, либо dt численного дифференцирования только переменной x2 и специального начального условия x2 (0) = x1 (0) для второй переменной.

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

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

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

3.2 Структура анализатора в пакете MVS.

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

3.2.1 Вычислимые последовательности формул.

Рассмотрим вектор s размерности n, а также вектор-функцию g ( s ) : n m, n m. Будем говорить об упорядоченных формулах или последовательности формул si = g i ( s );

i = 1, m, если формулы следует вычислять в заданном порядке. Разделим компоненты вектора s на два подвектора s1 и s 2, размерности n1, n2, n1 + n2 = n соответственно, и будем называть компоненты вектора s1 данными, а компоненты вектора s 2 вычисляемыми переменными. Назовем последовательность формул s k = g k ( s1, s 2 ), k = 1, n 2 вычислимой, если для каждого k в правую часть формулы входят только данные s1 или ранее вычисленные переменные s2.

вектора Наша задача определить, является ли конкретная последовательность формул вычислимой, а если нет, то можно ли переупорядочить ее так, чтобы она таковой стала.

Для ответа на эти вопросы введем матрицу структуры Q, столбцам которой соответствуют вычисляемые переменные, взятые в том же порядке, как они появляются в последовательности формул, а строкам – формулы последовательности. Пометим вхождение i-той компоненты вектора s 2 в k-тую формулу единицей в матрице структуры q k,i = 1, и нулем ее отсутствие - q k,i = 0.

Одним из возможных способов определение вычислимости последовательности формул связано с возможностью приведения ее к нижнему треугольному виду.

Пример 1. Вычислимая последовательность формул и ее матрица структуры Q.

2 s2 s s 2 1 s3 = f1 = s1 + s1 0 0 f s1 = f2 = s s 1 0 f s = f3 = s1 0 0 f Результат приведения матрицы структуры Q к нижнему треугольному виду.

2 s1 s s 2 1 s3 = f1 = s1 + s1 0 0 f s1 = f2 = s2 s 1 0 f 2 s1 = f 3 = s 0 1 f Вычислимая последовательность: f 3, f1, f Условия вычислимости последовательности формул можно сформулировать следующим образом. Для того чтобы последовательность формул была вычислимой (структурно невырожденной), должна существовать матриц перестановок P, такая, чтобы матрица P1 (Q) P1T имела бы блочно-диагональную структуру, и каждый блок являлся бы строго нижней треугольной матрицей (то есть матрицей с нулевой диагональю).

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

Конструктивный способ построить нужную форму может быть записан в виде алгоритма.

Обозначим через i номер строки матрицы структуры.

Алгоритм_3.1 определения вычислимости последовательности формул i = До тех пор, пока i n Найти все строки у которых q kj = 0, j = i, n и расположить их последовательно, начиная с i -той строки. При перемещении каждой найденной строки на новое место, а) размещенную там строку помещать на место найденной, б) менять местами столбцы с номерами, соответствующими старому положению найденной строки и новому.

Увеличить значение i на число перемещенных строк.

Конец цикла Этот же алгоритм может быть реализован с меньшими затратами.

Алгоритм_3.2 определения вычислимости последовательности формул i = До тех пор, пока i n Найти строку, у которой q kj = 0, j = 1, n i + 1, и поместить ее номер в массив p (i ) = k.

Вычеркнуть из матрицы k -тую строку и k -тый столбец, образовав новую матрицу меньшей размерности.

i = i + Конец цикла Завершение работы алгоритма приведет либо к построению блочной нижней треугольной, либо блочной матрицы. В первом случае будет построена вычисляемая последовательность формул, во втором – вычислимая последовательность максимально возможной длины и система уравнений.

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

3.2.1 Структурно невырожденные системы нелинейных алгебраических уравнений Рассмотрим вектор s размерности n и вектор-функцию v(s ) размерности m. Значение m не обязательно совпадает со значением n.

Будем называть систему v( s ) = 0 - формально определенной, если число уравнений равно числу неизвестных - n = m ;

переопределенной, если число уравнений больше числа уравнений - n m ;

и недоопределенной, n m.

если число уравнений меньше числа неизвестных Переопределенные и недоопределенные системы назовем структурно вырожденными. Рассмотрим задачу об определении структурной невырожденности формально определенной системы уравнений v( s ) = 0, где v вектор-функция размерности n.

Сопоставим формально определенной системе v( s ) = 0 квадратную матрицу структуры Q, строкам которой соответствуют уравнения системы, а столбцам – компоненты вектора s. Пометим вхождение i-той компоненты вектора s в k-тое уравнение единицей в матрице структуры q k,i = 1, и нулем ее отсутствие - q k,i = 0. Очевидно, что структурная невырожденность системы связана с рангом матрицы Q, и если rank (Q) = n, то система структурно не вырожденна.

Конструктивный способ определить структурную невырожденность системы заключается в возможности последовательно выражать компоненты вектора неизвестных si = g i ( s1,..si 1, si +1,...s n ), i = 1, n, начиная с первого уравнения системы, и подставлять в оставшиеся j = i + 1, n (прямой ход процесса подстановки) до тех пор, пока не будет построено уравнение, зависящее только от одной переменной. Для определения структурной невырожденности необязательно уметь явно строить явно определенную si = g i ( s1,..si 1, si +1,...s n ), подстановку нужна только информация о вхождении неизвестных в соответствующее уравнение.

Прямой ход процесса подстановки, можно описать последовательностью матриц Qi, i = 0, n 1;

Q0 = Q. Каждая из матриц, начиная с i = 1, совпадает по форме с матрицей, соответствующей исключению i-той неизвестной в методе Гаусса при решении систем линейных алгебраических уравнений. Единственным отличием процесса исключения в методе Гаусса и процесса проверки невырожденности является правило формирования новых строк в подматрице (i 1) (i 1) на каждом шаге.

При «подстановке» i-той неизвестной в оставшиеся j = i + 1, n уравнения при проверке невырожденности, в новом уравнении число единиц может только увеличится – к уже существующим единицам j-той строки добавляются не совпадающие по положению единицы i-той строки. Этот процесс в точности эквивалентен процессу появления новых ненулевых элементов в методе Гаусса при использовании любого ненулевого элемента в столбце в качестве ведущего, что позволяет использовать для определения невырожденности уже существующее программное обеспечение. Таким образом, если существует матрица перестановок P1, такая что матрицу P1Q нелинейной алгебраической системы описанным выше алгоритмом формальной подстановки можно привести к верхней треугольной матрице, то система уравнений структурно невырождена.

Пример 2. Система уравнений и ее матрица структуры, где - x, y, z, w неизвестные, ci - константы.

y x w z x + y c1 = 0 = f 1 1 1 0 f x y c2 = 0 = f 2 1 1 0 f x * y c3 = 0 = f 3 1 1 0 f z * w c4 = 0 = f 4 0 0 1 f Результат работы алгоритма псевдо-исключения показывает, что данная система структурно вырождена.

y x w z x + y c1 = 0 = f 1 1 1 0 f x y c2 = 0 = f 2 0 1 0 f x * y c3 = 0 = f 3 0 0 0 f z * w c4 = 0 = f 4 0 0 1 f Пример 3. Структурно-невырожденная система может быть численно вырожденной.

Исходная система и ее матрица структуры.

y x w z x + y c1 = 0 = f 1 1 1 0 f x + y c2 = 0 = f 2 1 1 0 f z * w c3 = 0 = f 3 0 0 1 f z * w c4 = 0 = f 4 0 0 1 f Результат псевдо-исключения, показывает что система и структурно невырождена, однако конкретный набор коэффициентов делает ее числено вырожденной.

y x w z x + y c1 = 0 = f 1 1 1 0 f x + y c2 = 0 = f 2 0 1 0 f z * w c3 = 0 = f 3 0 0 1 f z * w c4 = 0 = f 4 0 0 0 f 3.2.3 Системы нелинейных алгебраических уравнений с подстановками.

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

В языке MVL нет типа данных «уравнение» и «формула». Это объясняется следующим. Исходное представление пользователем выражений в виде уравнений и формул не всегда оказывается таким же в итоговой системе уравнений. Достаточно привести в качестве примера блок-схему кольца для двух блоков с формулами y = f (x), y = g (x), x1 = y 2, x2 = y1, чтобы убедиться, в том, что в итоге возникает система нелинейных алгебраических уравнений. В некоторых пакетах, все формулы превращаются в уравнения, но это не целесообразно. Если все-таки формула в результате оказывается уравнением, формируется предупреждение о том, что желание пользователя не может быть удовлетворено, и необходимо либо проверить описание задачи, либо согласится «с мнением» компилятора.

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

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

3.3.4 Системы дифференциальных уравнений.

ds, s, t ) = 0 преобразуется в пакете к форме Свободная форма F ( dt ds =z, dt F ( z, s, t ) = и одновременно находится структурный индекс полученной системы.

ds + B * s = f (t ) Уравнения с постоянными коэффициентами A * dt приводятся к форме ds1 ds + A12 2 + B11 s1 + B12 s 2 = f A, dt dt B21 s1 + B22 s 2 = f и одновременно находится индекс полученной системы.

ds + B (t ) * s = f (t ) Уравнения с переменными коэффициентами A(t ) * dt приводятся к форме ds1 ds + A12 (t ) 2 + B11 (t ) s1 + B12 (t ) s 2 = f A11 (t ), dt dt B21 (t ) s1 + B22 (t ) s 2 = f но находится только структурный ранг системы.

ds + g ( s ) = f (t ) приводятся к форме Уравнения A(t ) * dt ds1 ds + A12 (t ) 2 + g1 ( s1, s 2 ) = f A11 (t ) dt dt g 2 ( s1, s 2 ) = f и также находится только структурный ранг системы.

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

К основным проблемам поиска точки переключения следует отнести:

• Построение надежного алгоритма локализации очередной точки переключения с заданной точностью (глобальный поиск).

• Построение быстрого алгоритма поиска локализованной точки (локальный поиск).

Первая проблема связана с выбором величины шага модельного времени.

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

Смена поведения в MVS возможна после:

• истечения заданного времени, • получения сообщения, • наступления истинности заданного предиката.

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

Во входном языке пакета условие смены поведения можно записывать в виде:

• булевской функции Prb : boolean boolean • предиката Prg : arg of any MVS type boolean.

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

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

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

Рассмотрим в качестве примера три гибридных системы.

Пример А.

ds = a1 s + b1 ;

a1 = 1;

b1 = 1;

s (0) = ;

while y dt y = s ds = a 2 s + b2 ;

a1 = 1;

b2 = 1;

;

while y s (0) dt y = s Если этот пример рассматривать как пример гибридного автомата, имеющего периодическое решение, то в нем может только неправильно вычисляться только период.

Пример Б ds = a1 s + b1 ;

a1 = 1;


b1 = 1;

s (0) = ;

while y dt y = s ds = a 2 s + b2 ;

a1 = 1;

b2 = 1;

dt ;

else.

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

Пример В ds = a1 s + b1 ;

a1 = 1;

b1 = 1;

s (0) = ;

while y dt y = s ds = a 2 s + b2 ;

a1 = 1;

b2 = 1;

;

while y dt y = s И, наконец, в последнем примере хорошо видно, что точкой переключения является не интервал как было в предыдущих примерах, а действительно точка, а именно y = 0 и продолжение решения невозможно ни слева, ни справа.

В последнее время все чаще возникает необходимость записывать условие смены поведения в виде корня алгебраического уравнения, в общем случае нелинейного:

ds1 ds = f ( s1, t ) 2 = f ( s 2, t ).

dt dt when g ( s1 ) = Начнем рассмотрение блока поиска точки переключения с алгоритмов, предназначенных для локализации единственной точки переключения для уравнений с продолжаемым (за точку переключения) решением.

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

В данном случае в семействе Model Vision применяется метод деления отрезка пополам и в качестве нового значения модельного времени выбирается значение правой точки окончательного локализующего интервала. Новые начальные условия в этой точке вычисляются по значениям текущих переменных на левом конце окончательного локализующего интервала.

Такой подход вполне оправдан для случаев, когда следующим текущим состоянием является длительное состояние. Если же следующим текущим оказывается одно или несколько «мгновенных» состояний (предикат смены поведения истинен уже в начальной точке), то в гибридном времени, вместо отрезков нулевой длительности для непрерывного времени, соответствующим таким состояниям, появляются отрезки не нулевой длины, равной длине окончательного локализующего интервала.

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

Рассмотрим в качестве примера случай, когда смена поведения текущего s 0, состояния определяется сменой знака переменной s = (t t*), 0, а начальные условия нового текущего состояния выходной переменной y, y = (t t*),. В этом случае новые начальные условия могут сильно зависеть от ошибок в локализации точки переключения и значения выходной переменной (чувствительности уравнений для новых начальных условий по отношению к возмущению входных данных).

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

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

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

В этом случае, после того как точка переключения локализована, вместо неравенств, описывающих смену поведения, строятся уравнения, и применяется различные модификации метода Ньютона.

Уравнения строятся по следующим правилам:

- все логические выражения в форме неравенств для вещественных переменных и их отрицания приводятся к уравнениям. s ( x) 0 s( x) = 0.

- логические выражения, связанные логическим «И», рассматриваются как совместная система уравнений - логические выражения, связанные логическим «ИЛИ», рассматриваются как независимые системы уравнений. Уравнения сортируются по степени сложности (как мера сложности рассматривается порядок уравнений), и решаются последовательно, пока не будет найдено решение хотя бы одной системы.

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

3.3.3 Поиск в условиях, когда точка переключения определяется как корень уравнения на решении, описывающем поведение.

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

Рассмотрим особенности локализации корня в условиях непродолжаемости решения за точку переключения.

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

Одним из способов, предложенных в работе [10] применительно к поиску точки переключения, является переход от поиска корня нелинейного алгебраического уравнения, к поиску стационарного режима дифференциального уравнения, для которого положение равновесия и является корнем исходного алгебраического уравнения. В вычислительной математике [11] такой прием называется поиском корня методом установления. Очевидным достоинством этого приема является то, что автоматически обеспечивается невозможность перейти за точку переключения.

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

Пусть требуется найти корень функции g ( s ) = 0 на решении уравнения ds = f ( s ), s (0) = s 0. В написанной системе осуществим переход от dt независимой переменной t к новой независимой переменной h. Переход к новой независимой переменной и дифференцирование функции g (s) по независимой переменной приводит к системе:

ds dt = f ( s) ;

dh dh dg g s dt g dt = = f (s) ;

dh s t dh s dh dt Как видно из последней записи, зависимость пока произвольна, и dh ее можно выбрать, так чтобы корень исходного алгебраического уравнения оказался бы устойчивой стационарной точкой новой системы. Обозначим dt =R эту зависимость через R, положив dh ds = f ( s ) R, s ( 0) = s 0 ;

dh dg g s dt g = = f ( s ) R, g (0) = g ( s 0 );

dh s t dh s dt = R;

t (0) = dh Выбор функции R в виде R = ± g (s ) приводит к наиболее простой схеме ds = f ( s ) g ( s ), s (0) = s0 ;

dh g dg =± f ( s ) g ( s ), g (0) = g ( s0 );

s dh dt = g ( s );

t (0) = dh в которой, однако, необходимо дополнительно выбирать знак, для обеспечения устойчивости.

ds = as + b, s (0) = s 0 и Примером может служить исходная система dt предикат s c, c = const, которые превращаются в систему ds = ±(as + b)( x c), s (0) = s0 ;

dh dg = ± (as + b)( x c), g (0) = g ( s0 );

dh dt = ±( x c);

t (0) = dh g Выбор функции в виде R = g ( s ) обеспечивает устойчивость s [11] и приводит к уравнениям g ds = g ( s ) f ( s ), s (0) = s0 ;

s dh dg = g ( s ) f ( s ), g (0) = g ( s0 );

dh g dt = g ( s );

t (0) = s dh И, наконец, выбор функции в виде (для одной независимой g ( s ) g переменной) R =, = const приводит также к устойчивой f ( s ) s схеме из статьи [10]:

g ds = g ( s ), s (0) = s0 ;

s dh dg = g ( s ), g (0) = g ( s0 );

dh g g ( s ) dt = ;

t (0) = s f ( s ) dh 3.4 Построение графиков функций.

Автоматизация вычислений во многих программных средах стала настолько привычной, что мы подсознательно считаем результаты машинных вычислений безусловно правильными и начинаем их проверять в последнюю очередь. К числу интуитивно «безошибочных»

математических операций относится и такая «простая», как рисование графика гладкой функции.

В математических пакетах (Maple, Mathematica, Matlab, MathCad), как правило, существует отдельная операция построения графика гладкой функции – обычно называемая Plot, и операция одновременного решения уравнения и его визуализации (например, DEplot в Maple для решения дифференциального уравнения). С точки зрения пользователя, которому надо всего лишь увидеть и проанализировать поведение функции или решение уравнения, обе операции отличаются лишь способом задания исследуемой функции.

3.4.1 Проблемы, возникающие при построении графиков.

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

(t p ) В качестве тестовой функции выберем функцию f (t ) = t + A * e.

Напишем дифференциальное уравнение, которому она удовлетворяет, в (t p ) 2(t p ) dx = g (t ), g (t ) = f ' (t ) = 1, x(0) = f (0) форме A*e и будем dt рисовать графики функций f (t ) и x(t ) на промежутках [0, T ], T = 2, 4, A = 1;

p = 1;

= 0.01 с помощью различных при заданных значениях математических пакетов (Таблица 1).

Таблица 1. Тестируемые пакеты и процедуры пакет Построение Решение дифференциального графика уравнения и построение графика Mathematica Plot NDSolve Maple Plot Deplot Matlab Plot OdeSolver* Solver = ’45’;

‘23’… *) - В пакетах Maple и Mathematica численный метод, используемый для решения дифференциальных уравнений, является параметром процедуры.

По умолчанию в Maple используется RKF45, а в пакете Mathematica – Lsode. В пакете Matlab метод выбирается пользователем, c помощью суффикса Solver в имени процедуры, например Ode45. Процедура Ode использует тот же метод Рунге-Кутта-Фельдберга, что и в пакете Maple и рекомендуется как первоначальная для любой системы дифференциальных уравнений, свойства которой заранее неизвестны.

График функции f (t ) при значении T = 2 приведен на Рис.3. На Рис.


4 показаны одновременно графики функций f (t ) и g (t ).

1. 0. 0.7 0.8 0.9 1.1 1.2 1. Рис. 3 График функции f (t ) при значении T = 0.4 0.6 0.8 1.2 1.4 1. - - - Рис. 4. Графики функций f (t ) и ее производной g (t ).

Прекрасно зарекомендовавшие себя пакеты Maple, Mathematica и Matlab одинаково хорошо воспроизводят графики функций f (t ) и x(t ) при значении параметра T = 2. В тоже время, при воспроизведении графиков в пакете Mathematica, уже при T = 4, вместо ожидаемого «правильного»

графика (Рис. 5) 1 2 3 Рис. 5. «Правильный» график функции f (t ) при значении T = 4, построенный с помощью процедуры Plot пакета Mathematica мы наблюдаем только линейную составляющую функции x(t ) (Рис. 6) 1 2 3 Рис. 6. «Неправильный» график функции x(t ) при значении T = 4, построенный с помощью процедуры NDSolve пакета Mathematica Полученный результат заставляет, прежде всего, проверить, как изменится характер графиков функции x(t ), при смене численных методов. Результаты экспериментов приведены в Таблицах 2-4.

При значении T = 12 с поставленной задачей правильно справилась только процедура Ode23 пакета Matlab.

Таблица 2. Графики функции x(t ). Maple, T= метод график gear Правильный mgear Неправильный lsode Пик на графике есть, но он другой формы Rkf45 Правильный Dverk78 Правильный Таблица 3. Графики функции x(t ). Mathematica,. T= метод график Lsode Неправильный Adams Неправильный Gear Неправильный RungeKutta Неправильный Таблица 4. Графики функции x(t ). Matlab, T= метод график Ode45 Правильный Ode23 Правильный Ode113 Неправильный Ode15s Неправильный Ode23s Неправильный Ode23t Неправильный Ode23tb Правильный Исходя из того, что процедура Plot во всех пакетах успешно рисует T = 2, 4, 12, графики при любых значениях параметра можно предположить, что источником ошибок является численный метод, поставляющий исходную таблицу для аналогичной процедуры при решении уравнений. Действительно, например, для T = 4, алгоритм выбора шага интегрирования процедуры Ode113 пакета Matlab «не замечает» резкого роста производной в окрестности точки t = 1 и, при правильном вычислении каждой отдельной точки таблицы, неправильно воспроизводит поведение функции в целом (Таблица 5).

В пакете Maple, численные методы можно вызвать непосредственно для построения численного решения в заданной точке t = 1 и убедиться в некорректности именно алгоритма выбора шага при воспроизведении решения на больших промежутках. Как и следует ожидать, значения x(t ), при t = 1, оказываются вычисленными в этом случае с заданной по умолчанию точностью всеми имеющимися в пакете численными методами.

Таблица 5. Численное решение для функции x(t ), полученное процедурой Ode113 при T = x (t ) t 0.5181 0. 0.7181 0. 1.1181 1. 1.5181 1. 1.9181 1. Особое положение в этих экспериментах занимает пакет MathCad.

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

Однако винить только численные методы в неправильном построении графиков нельзя.

3.4.2 Аккуратные графики Рассмотрим гладкую вещественную функцию x(t ) одной вещественной переменной t на промежутке [0, T ]. Выберем сетку по оси ординат {~i : x(0) = ~1 ~2...~N = x(T )}. На практике ее задает либо x xx x пользователь, либо она выбирается автоматически, по оценкам максимального и минимального значения функции на рассматриваемом промежутке. Для аккуратного построения графика функции достаточно, по выбранной сетке, построить соответствующую сетку по оси абсцисс {t i : 0 = t1 t 2...t N = T }, и, тем самым, множество опорных точек {(t i, ~ (t i )), i = 1, N }, x таких, что i = 1, N 1 ~ (t i ) max( x(t )) & min( x(t )), t [t i, t i +1 ] ~ (t i +1 ). Число x x точек по оси абсцисс может оказаться недостаточным для восприятия рисунка, как графика гладкой функции, поэтому внутри любого интервала [t, t ] следует использовать дополнительные точки (t, ~ (t )), лишь бы x i + i k k значения ~ (t k ) не выходили за границы ~(t i ), ~(t i +1 ). Если цена вычисления x x x дополнительных точек высока, то вместо функции x(t ) можно использовать любую ее аппроксимацию, с меньшей стоимостью вычислений.

На Рис. 7 показана функция (пунктирная линия), для которой была правильно выбрана сетка по оси абсцисс, и ее график (сплошная линия) для заданных опорных точек. На Рис. 8 показана функция (сплошная линия), для которой неправильно выбрана сетка по оси абсцисс. Для построения ее графика (пунктирная линия), была использована интерполяционная функция. Как хорошо видно, в последнем случае мы сталкиваемся с тем же эффектом, что и при рисовании графиков решения дифференциального уравнения.

Алгоритм, рисующий график на заданной сетке, должен либо аккуратно строить график, либо сообщать о необходимости изменить сетку.

Рис. 7 Пример аккуратно нарисованного на заданной сетке графика Рис. 8. Пример неправильно построения опорных точек и неаккуратно нарисованного на заданной сетке графика.

Возможна более жесткая постановка задачи. По заданной сетке, алгоритм должен уметь выделять участки на оси абсцисс, где можно аккуратно построить график, и участки, где это сделать невозможно, и рисовать график так, как, например, делает это пакет Mathematica (Рис. 9).

7. 2. 0.25 0.5 0.75 1 1.25 1.5 1. -2. - -7. - Рис. 9. График производной функции f (t ), нарисованный процедурой Plot пакета Mathematica, при явно заданном ограничении на максимально допустимое значение по оси ординат.

3.4.3 Синхронное и асинхронное построение графиков Будем различать два способа построения графиков – синхронный (одновременный с построением таблицы) и асинхронный. При асинхронном способе, что характерно для математических пакетов, сначала строиться вся таблица необходимых опорных точек, а затем вызывается процедура построения графика. Обычно в этом случае при вычислении опорных точек используется сама функция x(t ), а дополнительные точки, необходимые для создания эффекта гладкости, находятся с помощью интерполирующей функции. При синхронном способе последовательно находится очередная точка сетки по оси абсцисс и немедленно прорисовывается очередной участок графика. Этот способ (его также иногда называют анимацией, хотя всего лишь рисуется след точки) характерен для пакетов, моделирующих динамические системы (Simulink, Model Vision Studium, AnyLogic). Последовательное вычисление точек по оси абсцисс в синхронном способе является очень важным, так как ось абсцисс в этом случае является осью времени.

Алгоритмы синхронного и асинхронного построения графиков могут различаться по виду используемой информации – алгоритмы, использующие только значения функции x(t ), алгоритмы, которым доступна также производная функции и так далее. Учитывая, что при моделировании динамических систем, функция x(t ) обычно является вектор-функцией (ограничимся вектором с двумя координатами x(t ) 2 ) и может быть задана различными способами, будем различать виды задания:

1. Явным способом в виде вектор-функции x(t ) = [ f1 (t ), f 2 (t )] 2. В виде системы дифференциальных уравнений dx dt = f1 ( x1, x2, t ), x1 (0) = x dx2 = f ( x, x, t ), x (0) = x dt 2 1 2 2 3. В виде системы алгебро-дифференциальных уравнений dx = f1 ( x1, x2, t ), x1 (0) = x dt 0 = f 2 ( x1, x2, t ), x2 (0) = x с согласованными (или несогласованными) начальными условиями по второй (алгебраической) компоненте.

4. В виде системы алгебраических уравнений 0 = f1 ( x1, x 2, t ), x1 (0) = x 0 = f 2 ( x1, x 2, t ), x 2 (0) = x с согласованными (или несогласованными) начальными условиями по обеим компонентам.

Помимо основной задачи построения графика функции x(t ), рассмотрим дополнительную задачу построения графика функции x1 = g ( x2 ), x1 = g ( x2 );

x1 = x1 (t );

x2 = x2 (t ), заданной параметрически:

необходимую для построения фазовых траекторий при исследовании поведения динамических систем.

То, что проблемы существуют и при синхронном рисовании графиков, подтверждает еще один эксперимент. В пакете Simulink «соберем»

функцию f (t ), рассмотренную во введении, из типовых блоков (Рис. 10) Рис. 10. Функция f (t ), реализованная средствами Simulink и воспроизведем ее для T = 2 (Рис. 11) и T = 4 (Рис. 12). И если раньше можно было грешить на процедуру автоматического выбора шага интегрирования, то теперь явно видна ошибка процедуры рисования графиков.

Рис. 11. Simulink. График функции f (t ),T= Рис. 12. Simulink. График функции f (t ), T= Аналогичная ситуация имеет место и при интегрировании дифференциального уравнения (Рис. 13):

Рис.13. Simulink. Дифференциального уравнения для функции x(t ) Решение x(t ) при T=4 даже методом Рунге-Кутта получается достаточно неточным (Рис. 14) Рис. 14 Асинхронное рисование графиков в пакете Simulink при решении дифференциального уравнения. Хорошо видна линейная интерполяция для промежуточных точек.

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

Литература 1. H. Elmqvist, S.E. Mattson, M. Otter. Modelica – a language for physical system modeling, visualizing and interaction. The 1999 IEEE Symposium on computer-aided control system design. CACSD,99. Hawaii, August 22-27, K. E. Brenan, S. L. Campbell, L. R. Petzold. Numerical solution of initial 2.

value problems in differential-algebraic equations. Elsever Science Publishers Co., Inc. 1989], [Uri. M. Ascher,Linda R. Petzold. Computer methods for ordinary differential equations and differential-algebraic equations. SIAM, P.I. Barton, The modeling and simulation of combined discrete/continues 3.

processes, PhD thesis, University of London, London, U. K., 1992) K. E. Brenan, S. L. Campbell, L. R. Petzold. Numerical solution of initial 4.

value problems in differential-algebraic equations. Elsever Science Publishers Co., Inc. А.Ф. Филиппов. Дифференциальные уравнения с разрывной правой 5.

частью, M.: Наука, 1985, 223 стр.

H. Elmqvist, S.E. Mattson, M. Otter. Modelica – a language for physical 6.

system modeling, visualizing and interaction. The 1999 IEEE Symposium on computer-aided control system design. CACSD’99. Hawaii, August 22-27, С. Писсанецки. Технология разреженных матриц. М.: Мир, 1988, 7.

стр.

C.W. Gear. Differential-algebraic equation index transformations.

8.

SIAM,J.Sci.Stat.Comp., 9(1988), 39- Constantios C. Pantelides. The consistent initialization of differential 9.

algebraic systems. SIAM J. Sci. Stat. Comput. Vol.9, № 2, March J. M. Esposito, V. Kumar, and G. Pappas. Accurate event detection for 10.

simulating hybrid systems.

Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы.

11.

М.: Наука, 1987,595 стр.

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

Первые версии численной библиотеки были созданы в те годы, когда создание DLL было трудоемкой операцией, поэтому оказалось проще переписать основные численные методы на язык Pascal. Библиотека пакета MV 2.1 была ориентирована на вариант гибридного автомата, который сегодня используется в пакете Stateflow.

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

Однако при создании пакета AnyLogic, написанного на языке Java, вновь стал вопрос о переводе исходных фортрановских текстов на новый язык.

Как показало выборочное тестирование, замедление работы при переходе на язык Java оказалось вполне приемлемым (не более чем в два, три раза), зато пакет теперь создавал «чистые» Java программы. Наличие прототипа библиотеки на языке FORTRAN позволило тщательно протестировать вновь созданную версию библиотеки. Вместе с библиотекой на язык Java были перенесены и программы тестирования, что позволило сравнивать сначала результаты работы FORTRAN библиотеки и JAVA библиотеки, а затем воспроизводить тесты одновременно в пакетах MVS и AnyLogic.

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

4.1 Организация библиотеки Схематически библиотеку численных методов MVS можно изобразить следующим образом (табл.1).

Таблица BLAS FLOAT LINPACK EISPACK MINPACK LA GNA NAE_SLV ODE_SLV DAE_SLV Программы библиотеки представляют собой иерархическую структуру. Уровни (строки таблицы) являются независимыми группами, использующими только программы предыдущих уровней (табл. 2).

Таблица BLAS Систематизированная коллекция подпрограмм [1], реализует основные элементарные операции линейной алгебры. Используется для структуризации и ускорения работы программ высшего уровня.

FLOAT Подпрограммы [2] этого раздела вычисляют основные характеристики машинной арифметики.

LINPACK Систематизированная коллекция подпрограмм [3] предназначена для решения систем линейных алгебраических уравнений с матрицами различного типа.

EISPACK Систематизированная коллекция подпрограмм [4] предназначена для решения проблемы собственных значений.

MINPACK Систематизированная коллекция подпрограмм [8] предназначена для нахождения минимумов функционалов и решения систем нелиненых алгебраических уравнений.

LA Подпрограммы этой группы реализуют основные матричные операции.

GNA Группа подпрограмм [5-6], предназначенная для решения основных задач численного анализа общего назначения.

NAE_SLV Группа подпрограмм, предназначенная для решения систем нелинейных алгебраических уравнений [7-8], возникающих при описании поведения.

ODE_SLV Группа подпрограмм, предназначенная для решения систем обыкновенных дифференциальных уравнений [9-10], возникающих при описании поведения.

DAE_SLV Группа подпрограмм, предназначенная для решения систем алгебро-дифференциальных уравнений [11], возникающих при описании поведения.

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

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

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

Исходя из сказанного, в MVS, пользователю предлагаются три категории программных реализаций численных методов.

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

Решатели для задач конкретного типа. Если пользователю известен тип решаемой задачи, то достаточно просто выбрать один из предлагаемых на выбор методов соответствующего класса. Однако и в этом случае, пользователь имеет дело не непосредственно с программной реализацией из соответствующей систематизированной коллекции, а с управляющей программой, анализирующей коды завершения вызываемой подпрограммы и, в случае необходимости, меняющей ее параметры так, чтобы обеспечить получение решения на всем интервале моделирования.

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

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

Аварийная остановка любого метода сопровождается диагностическим сообщением. В случае аварийного завершения работы, можно попытаться сменить метод, но лучше устранить причину прекращения счета. "Плохие" численные задачи, конечно же, существуют, но практика показывает, что плохо сформулированных задач гораздо больше.



Pages:     | 1 | 2 || 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.