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

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

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


Pages:   || 2 |
-- [ Страница 1 ] --

ФГБОУ ВПО «Воронежский государственный технический университет»

На правах рукописи

БЕЛЯНИН Алексей

Михайлович

РАЗРАБОТКА И РЕАЛИЗАЦИЯ ЭФФЕКТИВНЫХ ЧИСЛЕННЫХ

МЕТОДОВ МОДЕЛИРОВАНИЯ И ОПТИМИЗАЦИИ НА

ОСНОВЕ МЕТОДА МОМЕНТОВ

Специальность 05.13.18 – Математическое моделирование,

численные методы и комплексы программ

Д И С СЕ РТ АЦ И Я на соискание ученой степени кандидата технических наук

Научный руководитель — заслуженный деятель науки РФ, доктор технических наук, профессор Подвальный С.Л.

Воронеж 2014 Содержание Введение........................................................................................................................... Глава 1 Обзор эффективных численных методов моделирования и оптимизации на основе метода моментов.......................................................................................... 1.1 Метод моментов для построения математических моделей........................... 1.2 Методы решения дифференциальных уравнений............................................ 1.3 Численное решение задачи Коши...................................................................... 1.3.1 Одношаговые методы Рунге-Кутта................................................................. 1.3.2 Явные многошаговые методы......................................................................... 1.3.3 Неявные многошаговые методы..................................................................... 1.3.4 Формула дифференцирования назад.............................................................. 1.3.5 Сравнительный анализ методов решения обыкновенных дифференциальных уравнений................................................................................. 1.3.6 Оценка локальной погрешности решения...................................................... 1.4 Безградиентные методы оптимизации билинейных динамических систем.... 1.5 Градиентные методы оптимизации билинейных динамических систем....... 1.6 Метод модельных функций................................................................................ 1.7 Кинетически обоснованные модельные функции............................................ 1.8 Цель работы и задачи исследования.................................................................. Глава 2 Математическое моделирование билинейных динамических систем повышенной размерности на основе метода моментов............................................ 2.1 Численное сравнение модели кинетики в общем случае и с использованием метода моментов........................................................................................................ 2.2 Специфика моделирования билинейных динамических систем повышенной размерности................................................................................................................ 2.3 Сравнение методов решения систем обыкновенных дифференциальных уравнений.................................................................................................................... 2.4 Алгоритмы выбора длины шага................

......................................................... 2.5 Алгоритм переменного порядка и шага............................................................ 2.6 Алгоритм понижения порядка решаемой системы обыкновенных дифференциальных уравнений................................................................................. Выводы........................................................................................................................ Глава 3 Методы оптимизации билинейных динамических систем повышенной размерности на основе метода моментов................................................................... 3.1 Сравнение методов оптимизации....................................................................... 3.2 Алгоритм оптимизации билинейной динамической системы, состоящей из нескольких систем в моментах................................................................................. 3.2.1Алгоритм определения количества систем в моментах................................ 3.2.2 Выбор модельной функции............................................................................. 3.2.3 Определение начальных значений для решения........................................... 3.2.4 Определение математического ожидания для каждой системы в моментах..................................................................................................................... 3.2.5 Определение констант модели........................................................................ 3.2.6 Определение концентрации для каждой системы в моментах.................... 3.2.7 Решение обратной задачи................................................................................ 3.3 Численное решение задачи оптимизации билинейных динамических систем повышенной размерности на основе метода моментов......................................... 3.3.1 Пример билинейной динамической системы, состоящей из двух систем в моментах..................................................................................................................... 3.3.2 Пример билинейной динамической системы, состоящей из трх систем в моментах..................................................................................................................... Выводы........................................................................................................................ Глава 4 Разработка программного обеспечения для многоальтернативного моделирования и оптимизации билинейных динамических систем повышеной размерности.................................................................................................................... 4.1 Структура программного комплекса................................................................. 4.2 Модульная структура программного средства................................................. 4.3 Алгоритм работы программного модуля.......................................................... 4.3.1 Алгоритм определения количества систем в моментах............................. 4.3.2 Алгоритм определения начальных значений............................................... 4.3.3 Алгоритм решения задачи оптимизации...................................................... 4.4 Структуры базы данных.................................................................................... 4.5 Схема информационных потоков.................................................................... 4.6 Выбор среды разработки................................................................................... 4.7 Технические условия работы и запуск программы........................................ 4.8 Интерфейс программы...................................................................................... 4.9 Результаты работы программы........................................................................ Выводы...................................................................................................................... Основные результаты работы.................................................................................... Список использованной литературы......................................................................... ВВЕДЕНИЕ Актуальность темы.

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

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

Билинейные динамические системы повышенной размерности встречаются в процессах с цепными реакциями, такими как: процессы окисления (горение, взрыв), крекинга, полимеризации и другие. Цепные реакции применяются в химической и нефтяной промышленности. Моделирование таких процессов требует очень больших вычислительных затрат. Поэтому для решения билинейных динамических систем повышенной размерности используют метод моментов. Идея метода заключается в замене истинных соотношений выборочными аналогами. Метод моментов позволяет значительно снизить размерность решаемой билинейной динамической системы. Исследованием различных задач оптимального управления конечномерными системами с помощью метода L-проблемы моментов занимался Красовский Н.Н.

Теоретические вопросы моделирования динамических систем химической физики изучали Берлин А.А., Ениколопов Н.С. Исследованием процессов синтеза полимеров занимались Кафаров В.В., Подвальный С.Л., Спивак С.И., Будтов В.П.

Численные методы решения жестких систем изучали Новиков Е.А., Ракитский Ю.В., Черноруцкий И.Г.

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

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

Тематика диссертационной работы соответствует одному из научных направлений ФГБОУ ВПО «Воронежский государственный университет»

«Вычислительные комплексы и проблемно-ориентированные системы управления».

Цель работы и задачи исследования.

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

Для достижения поставленной цели необходимо решить следующие задачи:

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

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

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

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

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

Методы исследования.

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

Соответствие диссертации паспорту специальности.

Работа соответствует следующим пунктам паспорта специальности 05.13. «Математическое моделирование, численные методы и комплексы программ»:

П.3. Разработка, обоснование и тестирование эффективных вычислительных методов с применением современных компьютерных технологий.

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

П.8. Разработка систем компьютерного и имитационного моделирования.

Научная новизна.

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

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

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

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

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

Практическая значимость работы.

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

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

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

Реализация и внедрение результатов работы.

Основные алгоритмы и методы, предложенные в диссертации, реализованы и апробированы в виде программного комплекса моделирования и оптимизации билинейных динамических систем повышенной размерности на основе метода моментов. Система внедрена и используется в учебном процессе ФГБОУ ВПО «Воронежский государственный технический университет».

Апробация работы.

Основные положения диссертационной работы докладывались и обсуждались на Всероссийской научно-технической конференции «Перспективные исследования и разработки в области информационных технологий и связи» (Воронеж, 2012), Всероссийской конференции «Интеллектуальные информационные системы» (Воронеж, 2012), Международной молодежной научной школе «Теория и численные методы решения обратных и некорректных задач» (Воронеж, 2012), молоджной конференции «Интеллектуальные технологии будущего. Естественный и искусственный интеллект» (Воронеж, 2011), Всероссийской конференции «Новые технологии в научных исследованиях, проектировании, управлении, производстве» (Воронеж, 2013).

Публикации.

По результатам исследований опубликовано 8 научных работ, в том числе – в изданиях, рекомендованных ВАК РФ, получено 1 свидетельство на программу для электронных вычислительных машин, базу данных, топологию интегральных микросхем. В работах, опубликованных в соавторстве и приведенных в конце автореферата, лично автором предложены: [1, 2] – численные методы решения прямой и обратной кинетической задачи;

[3] — численные методы понижения порядка решаемой системы обыкновенных дифференциальных уравнений с переменным шагом;

[4, 5] — программный комплекс для создания специальных программных средств;

[6] — библиотека моделей на основе метода моментов;

[7, 8] — блок методов оптимизации при принятии решений.

Структура и объем работы.

Диссертационная работа изложена на 125 страницах, включает 27 таблиц и 39 рисунков;

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

Содержание работы.

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

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

Проанализированы явные и неявные, одношаговые и многошаговые методы решения обыкновенных дифференциальных уравнений (ОДУ): метод Рунге Кутта, метод Адамса-Башфорта, метод Адамса-Моултона, формула дифференцирования назад (ФДН). Выявлены трудности при решении систем ОДУ повышенной размерности.

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

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

экспоненциальное распределение, гамма-распределение, распределение Бизли.

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

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

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

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

Предложен оригинальный алгоритм выбора длины шага при решении систем обыкновенных дифференциальных уравнений с переменным шагом.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

В качестве интегрированной среды объектно-ориентированного программирования был выбран «Borland C++ Builder 6.0» с встроенным языком высокого уровня С++.

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

Разработан интерфейс программы.

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

ГЛАВА 1 ОБЗОР ЭФФЕКТИВНЫХ ЧИСЛЕННЫХ МЕТОДОВ МОДЕЛИРОВАНИЯ И ОПТИМИЗАЦИИ НА ОСНОВЕ МЕТОДА МОМЕНТОВ 1.1 Метод моментов для построения математических моделей В ряде случаев моделируемый процесс описывается билинейной динамической системой, состоящей из огромного количества обыкновенных дифференциальных уравнений (в общем виде из бесконечного числа уравнений).

Билинейные динамические системы повышенной размерности встречаются в процессах с цепными реакциями, такими как: процессы окисления (горение [9, 10], взрыв [11]), крекинга [12], полимеризации [13-16] и другие. Цепные реакции применяются в химической и нефтяной промышленности. Моделирование таких процессов требует очень больших вычислительных затрат. Поэтому для решения билинейных динамических систем повышенной размерности используют метод моментов. Идея метода заключается в замене истинных соотношений выборочными аналогами. Метод моментов позволяет значительно снизить размерность решаемой билинейной динамической системы.

В общем виде билинейные динамические системы имеют вид:

dxi ai xi bij xi x j (1.1) dt i1 i 1 j Рассмотрим метод на примере модели кинетики для полимеризации полимеров.

Для исследования процессов полимеризации используется модель кинетики процессов в виде системы обыкновенных дифференциальных уравнений большой (в общем случае бесконечной) размерности. С помощью такой универсальной модели можно получить модель любого типового кинетического модуля, выбирая соответствующие протеканию процесса значения констант скоростей элементарных стадий[17]. Например, ниже представлена универсальная модель полимеризации полимеров:

dS K d 3 S Pk dt k dS ' K d 3 S Pk K d 4 S ' M dt k dP Ki1P0 M dt dP M 1 ( Ki1P0 Ki 2 K p P Kt 2 P ( K d 1 K d 4 ) Pk K d 4 S ' ) 1 dt k Kt1P ( Kt 3 P Kt 4 P ) Pk K d 2 P M k K d 3 SP 1 1 1 1 k 1 k..........

..........

..........

..........

....

dPx K p M 1 ( Px 1 Px ) Px ( K t1 K t 2 M 1 ( K t 3 K t 4 ) Pk dt k K d 1 M 1 ) K d 3 SPx K d 2 ( M x Pk Px M k ), для 2 x k 1 k dM M 1 ( K i1 P0 K i 2 Pk ( K p K d 1 K t 2 K d 2 ) K d 4 S ' ) dt k P1 ( K t1 K d 2 M k K d 3 S K t 4 Pk ) k 1 k..........

..........

..........

..........

.....

x k dM x Px ( Kt1 K d 1M1 Kt 4 Pk K d 2 M k ) 0.5Kt 3 Pk Pxk dt k 1 k 1 k K d 3SPx Kt 2 M1Px1 K d 2 M x Pk, для 2 x k P0 — концентрация инициатора;

P1, P2, …, Px — концентрации активно растущих цепей длины 1, 2, …, x;

M1 — концентрация мономера;

M2, …, Mx — концентрации неактивных цепей длины 2, …, x;

S, S` — концентрации вспомогательных веществ;

t — время протекания процесса;

Ki1, Ki2 — константы скорости инициирования;

Kp — константа скорости роста цепи;

Kt1, Kt2, Kt3, Kt4 — константы скорости обрыва;

Kd1, Kd2, Kd3, Kd4 — константы скорости передачи цепи.

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

С увеличением размерности системы время, затраченное на ее решение, растет. В результате анализа графиков зависимости времени решения от количества уравнений можно сделать вывод, что это экспонента (рисунок 1.1). Из графика видно, что время, затраченное на решение системы 100 000 уравнений примерно равно 14 минутам. Следовательно, если экстраполировать график, то для решения системы из 2 000 000 уравнений (R=106) потребуется примерно часов. А при решении обратной кинетической задачи, приходится решать систему ОДУ около тысячи раз. Следовательно, время вычислений будет примерно месяцев, что неприемлемо.

Рисунок 1.1 – Зависимость времени решения системы ОДУ от количества уравнений в системе Статистическая теория подобных процессов исходит из предположения о возможности характеристики решения на основе моментов различных порядков [17]. Важной характеристикой решения являются соотношения дисперсии и математического ожидания. Для анализа вводятся понятия моментов, обычно применяемые в статистике и теории вероятностей для оценки распределения случайных величин. При этом, естественно, могут быть использованы различные моменты распределения - начальные и центральные, нормированные и ненормированные и т.д. Таким образом, порядок решаемой системы уменьшается до десятков уравнений.

1.2 Методы решения дифференциальных уравнений Для решения обыкновенных дифференциальных уравнений применяются аналитические, приближенные и численные методы [18].

Аналитические методы позволяют получать решения дифференциальных уравнений через элементарные или специальные функции в конечном виде и являются эффективным средством исследования уравнений, однако применимы лишь для ограниченного класса дифференциальных уравнений — линейных, с постоянными коэффициентами и др. Эти методы в основном реализованы в универсальных математических пакетах MathCad [19], Matlab [20] и других.

Однако в практических задачах они оказываются часто неприменимыми.

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

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

Метод конечных разностей основывается:

1) на замене непрерывной области определения решения D дискретным множеством точек, называемым сеткой h;

2) на замене непрерывных функций дискретными (сеточными), определенными на введенной сетке изменения аргумента;

3) на замене производных, входящих в уравнение, конечными разностями.

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

1.3 Численное решение задачи Коши Методы численного решения задачи Коши для обыкновенных дифференциальных уравнений dui ( x) f i ( x, u1, u2,...,un ), x x0, ui ( x0 ) ui( 0), i 1,2,...,n dx разделяются на два класса:

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

К этим методам относятся методы Рунге–Кутта и метод решения с помощью рядов Тейлора;

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

1.3.1 Одношаговые методы Рунге-Кутта Рассмотрим задачу Коши для ОДУ первого порядка dy f ( x, y ) dx y ( a ) y на интервале [a, b].

Разобьем этот интервал на m отрезков точками хi с постоянным шагом h.

Точное решение задачи Коши в точке xi+1 на любом из частичных интервалов [ xi, хi+1 ] можно представить в виде ряда Тейлора с центром в точке xi. Пусть функция f(x,y) имеет р +1 непрерывную производную по обоим аргументам. Тогда y ( p ) ( xi ) p y ' ' ( xi ) h O(h p 1).

y ( xi 1 ) y ( xi ) y ' ( xi )h h...

2! p!

Производные в данном выражении вычисляются по формулам:

y' ( x) f ( x, y) f ( x, y) f ( x, y) f ( x, y) f ( x, y) y' ' ( x) y' f ( x, y) x y x y f ( x, y) f ( x, y) y ' ' ' ( x) f ( x, y) ( x x y f ( x, y) f ( x, y) ( f ( x, y)) f ( x, y) y x y … ( p 1) ) ( y ( p 1) ) f ( x, y ).

y ( p ) ( x) (y x y Пренебрегая остаточным членом, получаем дискретное уравнение h 2 f ( xi, yi ) f ( xi, yi ) yi hf ( xi, yi ) [ f ( xi, yi )]...

yi x y h p p 1 f ( xi, yi ) p 1 f ( xi, yi ) p [ f ( xi, yi )] x p 1 y p p!

Введя обозначение h f ( xi, yi ) f ( xi, yi ) ( xi, y i, h ) f ( xi, y i ) [ f ( xi, yi )]...

x y h p1 p1 f ( xi, yi ) p1 f ( xi, yi ) p [ f ( xi, yi )] x p1 y p p!

перепишем это дискретное уравнение в виде yi 1 yi h( xi, yi, h).

Полученная формула определяет двухточечную явную разностную схему. По этой формуле можно последовательно определить все значения уi. Положив в формуле р = 1, получаем формулу явного метода Эйлера [21-24] yi 1 yi hf ( xi, yi ).

Итак, для получения разностных схем методом рядов Тейлора [25, 26] необходимо иметь аналитические выражения полных производных по х от функции. В том случае, когда получение полных производных достаточно громоздко или отсутствует аналитическое выражение для, применять ряды Тейлора затруднительно, поэтому используют методы Рунге Кутта [27, 28].

Идея методов Рунге-Кутта состоит в представлении дискретной задачи в виде yi 1 yi h ( xi, yi, h),0 i m 1, где функция ( xi, yi, h) была бы максимально близкой к ( xi, yi, h) и приближала бы отрезок ряда Тейлора с точностью, но не содержала бы производных от Функции.

Представим функцию ( xi, yi, h) в виде ( xi, y i, h) c1 f ( xi, y i ) c 2 f ( xi ha 2 y i b21 hf ( xi, y i )), где – неизвестные коэффициенты, которые определяются из условия близости функции и.

Разложим второе слагаемое в формуле в ряд Тейлора в окрестности h с центром в точке до членов порядка :

f ( xi, yi ) ( xi, yi, h) c1 f ( xi, yi ) c2 f ( xi, yi ) c2 ha x f ( xi, yi ) c2 hb21 f ( xi, yi ) O(h 2 ) y f ( xi, yi ) f ( xi, yi ) (c1 c2 ) f ( xi, yi ) c2 ha2 c2 hb21 f ( xi, yi ) O(h 2 ).

x yi Приравняв коэффициенты при одинаковых членах и выражениях и, получаем (c1 c2 ) 1;

c2 a2 0.5;

c2b21 0.5.

Следовательно c1 c2 0. a2 b21 1.

С учетом полученных коэффициентов запишем формулу Рунге-Кутта второго порядка:

.

Аналогичным способом можно получить формулы Рунге-Кутта более высоких порядков.

Порядок методов Рунге-Кутта р определяется числом вычислений функции f(x,y) на одном шаге.

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

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

1.3.2 Явные многошаговые методы Многошаговые методы появились гораздо раньше, чем методы Рунге-Кутта.

Впервые их получил Дж.К. Адамс еще в 1855 г. [30].

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

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

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

Среди многошаговых методов наибольшее распространение в практике вычислений получили методы Адамса [31, 32].

Общий подход к построению многошаговых методов состоит в следующем.

Пусть у(х)- точное решение задачи Коши.

.

Если это уравнение проинтегрировать на отрезке [ ], то получим, где.

Заменим подынтегральную функцию интерполяционным полиномом Лагранжа Р(х) с узлами интерполирования, :

.

Проинтегрировав этот полином в пределах от, до получаем формулу двухшагового метода h yi (3 f i f y ).

i 1 i Строя аналогичным образом полиномы более высоких порядков, получаем соответствующие многошаговые формулы разных порядков. Все формулы, полученные с использованием такого подхода, называются формулами Адамса Башфорта (табл. 1.1).

Таблица 1.1 – Формулы Адамса-Башфорта Порядок Формула метода y hfi y i 1 i h yi (3 f i f y ) 2 i 1 i h yi (23 f i 16 f 5 f y ) 3 i 1 i 1 i h y (55 f i 59 f 37 f 9 f y ) 4 i 1 i 24 i 1 i2 i h yi (1901 f i 2774 f 2616 f 1274 f y i 1 i 1 i2 i 1274 f ) i h yi (4277 f i 7923 f 9982 f 7298 f y i 1 i 1 i2 i 2877 f 475 f ) i i 1.3.3 Неявные многошаговые методы Важным классом методов решения жестких систем уравнений являются неявные многошаговые методы [33-37].

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

Формулы, полученные на основе рассмотренного подхода, носят название неявных формул Адамса-Моултона (табл. 1.2). Формула Адамса-Моултона первого порядка называется неявной формулой Эйлера. Для получения значения решения на очередном шаге при использовании неявных методов необходимо в общем случае решать нелинейное уравнение или систему нелинейных уравнений.

Таблица 1.2 – Формулы Адамса-Моултона Порядок Формула метода y hf y i 1 i i h yi ( f f) y 2 i 1 2 i 1 i h yi (5 f 8 fi f y ) 3 i 1 12 i 1 i h yi (9 f 19 f i 5 f f y ) 4 i 1 24 i 1 i 1 i h y 646 f i 264 f 106 f y (251 f i 1 i 720 i 1 i 1 i 19 f ) i h yi 1427 f i 798 f 482 f y (475 f i 1 i 1 i 1 i 173 f 27 f ) i 3 i 1.3.4 Формула дифференцирования назад Для вывода формулы дифференцирования назад предположим, что y'(хn+1) аппроксимируется линейной зависимостью 1k y / i y, то есть, n1 n1i h i y / ( y y y... y (1.2) ) n1 h 0 n1 1 n 2 n1 k n1k где k - порядок многочлена, аппроксимирующего y(х), т.е. y( x) x x 2... x k, h - шаг аппроксимации [38].

01 2 k В данном случае предполагается, что аппроксимация ведется с постоянным шагом.

Коэффициенты а0, a1, …,ак выбираются таким образом, чтобы формула (1.2) давала точное значение у'п+1 в точке хn+1 ( y n/ 1 y / ( xn1 ) ), когда у(х) – любой многочлен степени т(т=0,1, 2,...,k).

Выведем ФДН второго порядка. Учитывая (1.2), запишем y / ( y y y ), n1 h 0 n1 1 n 2 n а само решение y(х) представим в виде полинома второй степени x x x x y( x) n1 n1 (1.3) 0 1 2 h h Продифференцируем полученное выражение по х и вычислим значение полученной производной в точке хп+ x x / ( x) 1 2 xn 1 xn 2 h 2 h.

y h Вычислим также значения полинома (1.3) в точках хп+х, хп, xn-1 и подставим в выражение (1.2).

1 ( ( ) ( 2 4 )).

h 00 1 0 1 2 20 1 h После преобразований получаем:

1 (0 1 2 )0 (1 2 2 )1 (1 4 2 ) 2.

Неизвестные коэффициенты 0, 1, 2 находим, решая систему уравнений 0 1 2 2 1 4 2 0.

В результате получаем: 0 = -1,5;

1 = 2;

3 = -0,5.

Для метода порядка система уравнений для определения k-го коэффициентов 0, 1,…, k имеет вид:

0 1 1 1... 0 1 2...k 1 0 1 4...k 2 0 (1.4).................

......

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

Если дискретные точки распределены на отрезке интегрирования неравномерно, то коэффициенты системы уравнений (1.4) уже будут функциями дискретных моментов времени, поэтому при переходе к новым точкам хп+2, хп+3...каждый раз нужно решать систему (1.4) заново.

Итак, для нахождения решения ОДУ на п+1 шаге с использованием ФДН нужно решить систему нелинейных уравнений вида y/, x ) 0, F( y n1, n1 n где y / ( x ) ( y y y... y ) n1 h 0 n1 1 n 2 n1 k n1k.

При использовании формулы дифференцирования назад время решения можно уменьшить, если для k вычисленных точек построить полином степени k, проходящий через k+1 точку (хп, хn-1,...,хп-к ) и с его помощью спрогнозировать первое приближение к yn+1. Обобщенное выражение для спрогнозированного k значения имеет вид y ( p) i y.

n1 i1 n1i Рассмотрим случай, когда k=2, то есть, y ( p) y n y y. (1.5) n1 1 2 n1 3 n Выберем полином второй степени в виде x x x x y( x) n1 n1 01 h h и найдем его значение в точках хn+1, xn, xn-1, xn-2.

y( x ), n1 y( xn ), y( x ) 2 4, n1 0 1 ) 3 9.

y( x n2 0 1 Полученные значения подставляем в (1.5) и после преобразований получаем 0 (1 2 3 )0 (1 2 2 3 3 )1 (1 4 2 9 3 ).

Коэффициенты находят из решения системы уравнений:

1 1 2 2 3 3 1 4 2 9 3.

Для случая произвольного k коэффициенты i определяются из системы вида:

0 1 1 1... 0 1 2...k 1 0 1 4...k 2....................

... 0 1 2 k...k k k Итак, общий алгоритмом решения любой системы неявных дифференциально-алгебраических уравнений с использованием ФДН можно представить так:

Шаг 1. Задать xn+1 =xn+h.

Шаг 2. Определить коэффициенты i (i 1,2,..., k 1) при х = хп+1 и вычислить спрогнозированное значение по формуле k y ( p) i y.

n1 i1 n1i Шаг 3. Определить коэффициенты i (i = 0,1,..., k) при х = хn+1 и решать, y/, x ) уравнение начальным приближением, равным F( y c n1 n1 n k 1 h y ( p) i y yp.

, методом Ньютона-Рафсона y n1 i1 n1i n1 n x x nk n1 Шаг 4. Вычислить ошибку.

Шаг 5. Сравнить E с заданным пользователем уровнем ошибки. Выбрать новый шаг по времени h и новый порядок k таким образом, чтобы значение h было максимальным из всех, для которых Е сравнимо с указанной пользователем максимальной ошибкой. Если Е очень велико, то отказаться от последнего шага.

Шаг 6. Выбрать п = п + 1 и перейти к шагу 1.

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

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

1. Явный одношаговый метод Эйлера.

2. Явный одношаговый метод Рунге-Кутта 4-ый порядок.

3. Явный многошаговый метод Адамс-Башфорт 3-ий порядок.

4. Неявный одношаговый метод Эйлера.

5. Явный многошаговый метод Адамса-Моултона 3-ий порядок.

6. Жестко устойчивый метод с использованием формулы дифференцирования назад 4-ый порядок.

В таблице 1.3 указаны основные характеристики для сравнения методов и даны номера методов, лучших с точки зрения этих характеристик.

Таблица 1.3 – Сравнение методов решения ОДУ Нежесткие системы Жесткие системы ОДУ ОДУ Устойчивость Не важна 6,5, Точность 2,3 6, Эффективность 2,3 Простота 1 Отсутствие сложностей при изменении 1,2 текущего шага Следует отметить эффективность многошаговых методов по сравнению с одношаговыми по затратам процессорного времени, но и их важный недостаток – необходимость выполнения первых нескольких шагов каким-либо одношаговым методом («разгон»). Оценка «не важна» в первой строке таблицы означает, что для нежестких систем при правильно выбранных шагах интегрирования устойчивость не нарушается.

1.3.6 Оценка локальной погрешности решения Для оценки локальной погрешности решения при переходе от точки xn к точке xn+1 используем также способ Рунге. В этом случае по одной и той же формуле вычисляются два приближения к решению в одной точке, но с разными шагами, которые используются для получения оценки.

Для локальной (в соседней точке) ошибки решения с шагом h справедлива оценка u ( xn1 ) y n1 ( xn, y n )h p 1 O (h p 2 ) При этом считается n = 0 — начальная погрешность.

Для локальной ошибки решения с шагом h/2 при xn+1/2 это выражение принимает вид h u ( xn h / 2) ~n1/ 2 ~ ( xn, yn )( ) p y При выполнении второго полушага из точки xn+h/2 в точку xn+1 получаем приближенное значение в точке хn+1 = xn+h:

h u ( xn h) ~n1 ~ ( xn h / 2, yn1/ 2 )( ) p1.

y Так как при расчетах используются малые величины h, то, в силу предположения линейности приращения функции, погрешность решения в точке хn+1 = xn+h (после второго полушага) приближенно можно оценить посредством выражения h u ( xn h) ~n1 ~ 2 ( xn, yn )( ) p1.

y На основании последних двух выражений получаем:

~ p 1 ~ ( y n 1 y n 1 ) ( xn, y n )h (1 1 / 2 p ) h p 1 ~ ( ~n 1 y n 1 ) y 2 ( x n, y n )( ).

( 2 p 1) Если в качестве приближения к решению в точке хn+1 принять, то погрешность метода равна ~ ~ u ( x h) ~ ~ ( y n 1 y n 1 ).

y n p n (1 1 / 2 p ) n Если же приближением считать, полученное с шагом h/2, то погрешность метода равна ~ ~ u ( x h) ~ ~ ( y n 1 y n 1 ).

y n p n (2 p 1) n 1.4 Безградиентные методы оптимизации билинейных динамических систем Безградиентные методы или методы нулевого порядка не требуют знания целевой функции в явном виде [39-47]. Они не требуют регулярности и непрерывности целевой функции и существования производных. Это является существенным достоинством при решении сложных технических и экономических задач [48]. При реализации прямых методов существенно сокращается этап подготовки решения задачи, так как нет необходимости в определении первых и вторых производных. К прямым методам относится целый ряд алгоритмов, которые отличаются по своей эффективности. Такие методы носят в основном эвристический характер. Прямые методы предназначены для решения безусловных задач оптимизации вида min f ( x).

xE n Рассмотрим безградиентный метод Гаусса. Это простейший алгоритм [49], заключающийся в том, что на каждом шаге (каждой итерации) минимизация осуществляется только по одной компоненте вектора переменных x.

Пусть нам дано начальное приближение x ( x10, x2,..., xn )T. На первой 0 итерации находим значение минимума функции при изменяющейся первой координате и фиксированных остальных компонентах, т.е.

x11 arg min f ( x1, x2,..., xn ).

0 x В результате получаем новую точку x ( x1, x2,..., xn )T. Далее из точки x 1 0 ищем минимум функции, изменяя вторую координату и считая фиксированными все остальные координаты. В результате получаем значение x1 arg min f ( x11, x2, x30,..., xn ) x и новую точку x ( x1, x2, x30,..., xn )T. Продолжая процесс, после n шагов получаем 1 2 n точку x ( x1, x2,..., xn )T, начиная с которой процесс поиска возобновляется по 1 2 n первой переменной.

В качестве условий прекращения поиска можно использовать следующие два критерия:

k ) f (x ) k 1) f ( x 2) xik 1 xik i, i.

Метод очень прост, но не очень эффективен. Проблемы могут возникнуть, когда линии уровня сильно вытянуты и «эллипсоиды» ориентированы, например, вдоль прямых вида x1 x2. В подобной ситуации поиск быстро застревает на дне такого оврага, а если начальное приближение оказывается на оси «эллипсоида», то процесс так и останется в этой точке [50-52].

Хорошие результаты получаются в тех случаях, когда целевая функция представляет собой выпуклую сепарабельную функцию вида n f ( x) f i ( xi ).

i Рисунок 1.2 – Пример траектории спуска в алгоритме Гаусса 1.5 Градиентные методы оптимизации билинейных динамических систем К градиентным методам или методам первого порядка относятся алгоритмы, в которых в процессе поиска кроме информации о самой функции используется информация о производных первого порядка. К группе таких методов относятся различные градиентные методы.

Рассмотрим простейший градиентный метод – метод наискорейшего спуска.

Градиент функции в любой точке показывает направление наибольшего локального увеличения f (x). Поэтому при поиске минимума f (x), следует двигаться в направлении противоположном направлению градиента в данной точке, то есть в направлении наискорейшего спуска [53].

Итерационная формула процесса наискорейшего спуска имеет вид k x k f ( x ), k k x или k f ( x ) k x x k S.

k k k k x k f ( x ) Очевидно, что в зависимости от выбора параметра траектории спуска будут существенно различаться (рисунок 1.3). При большом значении траектория спуска будет представлять собой колебательный процесс, а при слишком больших процесс может расходиться. При выборе малых траектория спуска будет плавной, но и процесс будет сходиться очень медленно.

Рисунок 1.3 – Траектории спуска в зависимости от величины шага Обычно выбирают из условия k arg min f ( x S ), k k решая одномерную задачу минимизации с использованием некоторого метода. В этом случае получаем алгоритм наискорейшего спуска.

Если определяется в результате одномерной минимизации, то градиент в точке очередного приближения будет ортогонален направлению предыдущего k k спуска f ( x ) S.

Вообще говоря, процедура наискорейшего спуска может закончиться в стационарной точке любого типа, в которой f ( x) 0. Поэтому следует проверять, не завершился ли алгоритм в седловой точке.

Эффективность алгоритма зависит от вида минимизируемой функции.

Алгоритм наискорейшего спуска сойдется за одну итерацию при любом начальном приближении для функции f ( x) x12 x2 (рисунок 1.4). Но сходимость будет очень медленной, например, в случае функции вида f ( x) x12 100x2. В тех ситуациях, когда линия уровня минимизируемой функции представляет собой прямолинейный или, хуже того, криволинейный «овраг» эффективность алгоритма оказывается очень низкой.

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

В алгоритме наискорейшего спуска на каждом этапе поиска используется k k только текущая информация о функции f ( x ) и градиенте f ( x ). В алгоритмах сопряженных градиентов используется информация о поиске на предыдущих этапах спуска.

Рассмотрим метод сопряженных градиентов.

k Направление поиска S на текущем шаге k строится как линейная k комбинация наискорейшего спуска f ( x ) на данном шаге и направлений k 0 спуска S, S,..., S на предыдущих шагах. Веса в линейной комбинации выбираются таким образом, чтобы сделать эти направления сопряженными. В этом случае квадратичная функция будет минимизироваться за шагов n алгоритма.

При выборе весов используется только текущий градиент и градиент в предыдущей точке.

0 В начальной точке x направление спуска S f ( x ) :

x x 0 S, 1 0 (1.6) где 0 выбирается из соображений минимальности целевой функции в данном направлении 0 arg min f ( x S ).

0 Новое направление поиска S f ( x ) 1 S, 1 1 (1.7) где 1 выбирается так, чтобы сделать направления S и S сопряженными по 1 отношению к матрице H :

0T S H S 0. (1.8) Для квадратичной функции справедливы соотношения:

1T 0 0 f ( x ) f ( x ) T f ( x ) x x 2 f ( x ) x, где x x x, 0 f ( x) f ( x ) 2 f ( x ) x.

Если положить x x, тогда x x 0 S и 1 1 0 f ( x ) f ( x ) 2 f ( x ) x x H0 S.

1 0 0 1 0 (1.9) 0T Воспользуемся (1.9), чтобы исключить S из уравнения (1.8). Для квадратичной функции H H T, поэтому транспонировав (1.9) и, умножив справа на H 1, получим 0T 1 0 [f ( x ) f ( x )]T H 1 x x H T H и далее 1 0 1 ( x x )T [f ( x ) f ( x )]T H 0T S 0 0.


Следовательно, для сопряженности S и S [f ( x ) f ( x )]T H 1H S [f ( x ) f ( x )]T [f ( x ) 1 S ] 0.

1 0 1 1 0 1 Вследствие изложенных ранее свойств сопряженности все перекрестные 0 члены исчезают [54-58]. Учитывая, что S f ( x ) и, следовательно, T f ( x )f ( x ) 1T f ( x )f ( x ) 0, 1 1 0 получим для 1 следующее соотношение:

f ( x ) 1 f ( x )f ( x ) T 1. (1.10) 0 0 T f ( x )f ( x ) f ( x ) Направление поиска S строится в виде линейной комбинации векторов 2 0 1 f ( x ), S, S, причем так, чтобы оно было сопряженным с S.

2 Если распространить сделанные выкладки на S, S,…, опуская их k kT k k содержание и учитывая, что S f ( x ) 0 приводит к T f ( x )f ( x ) 0, можно получить общее выражение для k :

k f ( x ) k. (1.11) k f ( x ) Все весовые коэффициенты, предшествующие k, что особенно интересно, оказываются нулевыми.

Полностью алгоритм описывается следующей последовательностью действий:

0 1. В точке начального приближения вычисляется S f ( x ).

2. На k-м шаге с помощью одномерного поиска в направлении определяется минимум функции, то есть решается задача k arg min f ( x S ) k k k x k S.

k k и находится очередное приближение x k 1 k 3. Вычисляется f ( x ) и f ( x ).

k 1 k f ( x ) k 1 S.

k 4. Определяется направление S 5. Алгоритм заканчивается, если S, где – заданная величина.

k После n 1 итераций ( k n ), если не произошел останов алгоритма, n процедура циклически повторяется с заменой x на x и возвратом на первый пункт алгоритма. Если исходная функция является квадратичной, то ( n 1)-ое приближение даст точку экстремума данной функции. Описанный алгоритм с построением k по формулам (1.10) соответствует методу сопряженных градиентов Флетчера-Ривса [59].

Рисунок 1.5 – Траектория спуска метода сопряженных градиентов для квадратичной функции В модификации Полака-Рибьера (Пшеничного) метод сопряженных градиентов отличается только вычислением k k k T f ( x )[f ( x ) f ( x )] k k 1 k T f ( x ) S.

В случае квадратичных функций обе модификации примерно эквивалентны.

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

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

Арсенал методов, используемых исследователями для этих целей, достаточно широк. Если известна кинетическая схема процесса, описывающаяся системой дифференциальных уравнений, то в принципе, для того, чтобы построить графическую кривую решения, нет необходимости знать аналитическое выражение для функции распределения. Для этой цели придется провести численное интегрирование большого количества [60-63] дифференциальных уравнений, которые выводятся исходя из условий материального баланса по каждому из компонентов реакции с учетом закона действующих масс. Решить подобную систему, не вводя никаких упрощающих предположений, не всегда представляется возможным. Для упрощения таких систем уравнений большое распространение получили метод «квазистационарного» состояния [64-67] и метод моментов [68-70].

Ввиду вероятностной природы подобных процессов, были развиты статистические методы анализа кривых решения, среди которых хотелось бы отметить аксиоматический подход, созданный С. Я. Френкелем [71] и дополненный и обобщенный в работах Шаманина [72]. Согласно этому подходу исследование каждой системы начинается с построения детальной вероятностной модели изучаемого процесса. При этом принимается факт общности экспоненциального распределения, сформулированного в виде отдельной теоремы. Причина подобной общности легко объяснима, если учесть, что всякий происходящий в природе процесс характеризуется неким параметром со средним временем релаксации, которое, как правило, имеет экспоненциальное распределение соответствующей компоненты.

Для решения обратной задачи формирования решения в технологических процессах различной природы наиболее интересен метод модельных функций [73]. Модельные функции применяются для аналитической аппроксимации экспериментальных кривых решения. Подбор модельной функции для описания экспериментально полученной кривой решения проводится следующим образом:

измеряются или вычисляются из экспериментальной кривой распределения дисперсия и математическое ожидание;

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

вычисляются параметры этих функций;

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

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

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

Большинство из них получено эмпирическим путем с целью удовлетворительного описания экспериментальных результатов. Но некоторые из модельных функций имеют четкий кинетический смысл. К ним относятся распределения Флори [74 76], Шульца-Флори [77], Пуассона [78] и Бизли [79]. Они были получены либо непосредственно при рассмотрении кинетики процесса, либо при анализе кинетической схемы с использованием статистического подхода [80-82].

В таблице 1.4 приведены характеристики наиболее известных модельных функций. Распределение Флори и Пуасона зависит только от одного параметра и требует наименьших вычислительных затрат, однако область применения у этих модельных функций сильно ограничена. Распределение Шульца-Флори и Бизли имеет большую область применения, однако зависит от двух параметров и требует больших вычислительных затрат.

Таблица 1.4 Характеристики различных модельных функций Распреде Распределение Распределение Распределение ление Шульца-Флори Пуасона Бизли Флори (M ) pe Дифференци M e ( ) (M ) альное (M ) e (1 p)2 Me численное M!

(1 M ) распределение 0, при M0, k i e Интегральное 1 [1 M (1 p)] i 0 i!, 1 e M численное e M распределение при (1 M ) k M k 1, k 0,1,2...

Отношение 2 p моментов первого и (1 ) нулевого порядков Отношение 64p моментов 1+ второго и (2 p ) (1 2 ) первого порядков 2 p 2, при 0p 1 Дисперсия 2 (1 ) 2 (1 2 ) 2 2, при p= 2, при p= 2 p 2 p, при 0p Среднечисле- нное 1 1, при p=1 отклонение 2, при p= 1.8 Цель работы и задачи исследования Целью исследования является разработка математических моделей, алгоритмов и эффективных численных методов для моделирования и оптимизации билинейных динамических систем повышенной размерности на основе метода моментов.

Для достижения поставленной цели необходимо решить следующие задачи:

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

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

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

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

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

ГЛАВА 2 МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ БИЛИНЕЙНЫХ ДИНАМИЧЕСКИХ СИСТЕМ ПОВЫШЕННОЙ РАЗМЕРНОСТИ НА ОСНОВЕ МЕТОДА МОМЕНТОВ Знание кинетических уравнений реакции в дифференциальной форме позволяет определить время достижения некоторой заданной концентрации реагирующего вещества (или продукта реакции). К сожалению, кинетическое уравнение реакции может быть получено только при е экспериментальном изучении. Поэтому важную роль для решения кинетических задач играет ЭВМ [1].

Для моделирования билинейных динамических систем повышенной размерности на основе метода моментов, в зависимости от степени ее жесткости, используются как явные, так и неявные численные методы решения систем нелинейных ОДУ [83]. Для работы с жесткими система ОДУ используют неявные многошаговые методы, например метод Адамса-Моултона или ФДН [38].

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

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

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

dxi ai xi bij xi x j.

dt i 1 i 1 j Рассмотрим следующую билинейную динамическую систему большой (в общем случае бесконечной) размерности для цепного процесса:

dP K i M 1 K p M 1 P1 K t P dt....................

........................

dPx K p M 1 ( Px 1 Px ) K t Px, 2 x dt dM K i M 1 K p M 1 Pk K t P dt k....................

.........................

(2.1) dM x K t Px dt где P1, P2, …, Px — активно растущие цепи длины 1, 2, …, x;


M1, M2, …, Mx — неактивно растущие цепи длины 2, …, x;

Ki, Kp, Kt — параметры модели.

Обычно решение представляется на фазовой плоскости qw*M для разных моментов времени. На рисунке 2.1 представлена зависимость кривой решения от количества уравнений в билинейной динамической системе (2.1) при Ki= 0,00005;

Kt=0,1;

Kp=32. Как видно из графика и таблицы 2.1, чем больше уравнений в билинейной динамической системе (2.1), тем более точная кривая решения. В данном случае значение математического ожидания равно 196,15. В реальных задачах это значение порядка 103–105. Соответственно, чтобы просчитать это 104– значение точно, необходимо решить уравнений. Проведнный эксперимент показывает, что время для решения общей системы уравнений кинетики процесса будет составлять примерно 5 часов, а для решения обратной задачи 7 месяцев, что неприемлемо.

Рисунок 2.1 – Кривая решения в зависимости от количества уравнений в билинейной динамической системе. Здесь математическое ожидание равно 196,15, число уравнений до Воспользуемся методом моментов. Моменты для билинейной динамической системы (2.1) вычисляются по следующим формулам:

j x j Px ;

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

dP K i M 1 ( K t K p M 1 ) P dt dM K i M 1 K p M 1 0 ( K t K p M 1 ) P dt d K p M 1 P1 K t dt d K p M 1 (2 P1 m 0 ) K t dt d K p M 1 (4 P1 2 1 0 ) K t dt d Kt dt d K t dt d Kt dt Построим решение, полученное с помощью метода моментов и модельной функции Флори ( qW ( M ) 2 Me M ), и сравним с полученными выше данными (Рисунок 2.2). Как видно из таблицы 2.2 использование метода моментов позволяет получать достаточно точное решение, при этом количество уравнений уменьшилось с нескольких тысяч до 8.

Рисунок 2.2 – Сравнение кривых решения, полученных решением системы в моментах и в общем виде Таблица 2.1 – Зависимость средних молекулярных масс от количества уравнений в модели кинетики Отношение Отношение Математи- моментов Математи- моментов Количество Количество ческое второго и ческое второго и уравнений уравнений ожидание первого ожидание первого порядков порядков 1100 160,3309 271,1709 2700 195,9073 386, 1200 166,3792 286,5639 2800 196,1404 388, 1300 171,5696 300,5648 2900 196,3171 389, 1400 175,995 313,2177 3000 196,4487 389, 1500 179,7511 324,5834 3100 196,5478 390, 1600 182,9153 334,7178 3200 196,6185 391, 1700 185,5691 343,6944 3300 196,6733 391, 1800 187,7759 351,5802 3400 196,7104 391, 1900 189,5998 358,4541 3500 196,7367 391, 2000 191,0973 364,3976 3600 196,7583 392, 2100 192,317 369,4915 3700 196,7708 392, 2200 193,3026 373,8185 3800 196,779 392, 2300 194,0932 377,4611 3900 196,7839 392, 2400 194,7233 380,5 4000 196,7858 392, 2500 195,2202 383,0097 4100 196,7862 392, 2600 195,6072 385,0602 4200 196,7862 392, Таблица 2.2 – Сравнение полученных результатов решения методом моментов и исходной системы Относительная ошибка отношение моментов Отношение времени вычисления Количество второго и первого исходной системы к времени уравнений порядков при решении вычисления методом моментов исходной системы, % 1100 18,5 1500 8,7 2000 2,9 2500 0,8 3000 0,2 3500 0,025 4000 0,0002 Специфика моделирования билинейных динамических систем 2. повышенной размерности Для билинейных динамических систем, состоящих из одной системы в моментах, модельные функции либо теоретически обоснованы, либо получены экспериментально. К теоретически обоснованным модельным функциям относят:

( M ) 1. Экспоненциальное распределение (Флори) qn ( M ) e.

( M ) ( M ) (1 p ) 2 Me 2. Распределение Шульца-Флори: qn ( M ) p e.

3. Распределение Бизли: qn ( M ).

(1 M ) Теоретически вид модельных функций доказан только для простейших моделей. Требуется численное обоснование вида модельной функции прямыми расчетами.

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

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

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

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

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

Сравним следующие методы для решения систем обыкновенных дифференциальных уравнений:

1. Явный одношаговый метод Рунге-Кутта 4-ый порядок.

2. Явный многошаговый метод Адамс-Башфорт 3-ий порядок.

3. Неявный многошаговый метод Адамс-Моултон 3-ий порядок.

4. Жестко устойчивый метод с использованием ФДН 4-ый порядок.

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

Модель №1 – случайный обрыв [69];

Модель №2 – мономерный обрыв;

Модель №3 – обрыв диспропорционированием;

Модель №4 – полимеризация бутадиена на каталитической системе NdCl3*3ТБФ-ТИБА [84,85];

Модель №5 – случайный обрыв с изменяющимися во времени параметрами;

Модель №6 – мономерный обрывом с изменяющимися во времени параметрами;

Модель №7 – обрыв диспропорционированием с изменяющимися во времени параметрами;

Для решения моделей был разработан программный комплекс, написанный на объектно-ориентированном языке программирования С++.

На отрезке tє [0, 10] все четыре метода при шаге h=0,00001 дают примерно одинаковые значения. В качестве эталона возьмм решение методом ФДН четвртого порядка при шаге h=0,00001.

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

1. Берм по 100 контрольных точек с эталона и решения исследуемого метода.

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

.

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

, где n – порядок химической системы, k – количество контрольных точек.

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

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

Ниже приведено сравнение методов решения систем обыкновенных дифференциальных уравнений по устойчивости и при tє[0, 10] (Таблица 2.3 - 2.6).

Из полученных данных можно сделать вывод, что наиболее устойчивым методом является метод ФДН.

Таблица 2.3 – Максимальный шаг, при котором критерий точности не превышает Модель Модель Модель Модель Модель Модель Модель №1 №2 №3 №4 №5 №6 № Рунге- 0,01923 0,01315 0,00316 0,00584 0,00769 0, 0, Кутт 1 8 5 8 2 Адамс 0,00259 0,00112 0,00645 0, Башфор 0,002398 0,008 1 4 2 т Адамс- 0,00892 0,01190 0,00724 0,00090 0, 0,012821 0, Моултон 9 5 6 9 0,09090 0,03703 0, ФДН 0,008475 0,125 0,125 0, 9 7 Таблица 2.4 – Максимальный шаг, при котором критерий точности не превышает Модель Модель Модель Модель Модель Модель Модель №1 №2 №3 №4 №5 №6 № Рунге 0,012048 0,111111 0,013514 0,003086 0,005882 0,009346 0, Кутт Адамс 0,002415 0,010309 0,002625 0,000667 0,001142 0,009709 0, Башфорт Адамс 0,016949 0,03125 0,013699 0,012195 0,008403 0,001468 0, Моултон ФДН 0,012195 0,125 0,1 0,25 0,071429 0,021277 0, Таблица 2.5 – Ранжирование методов при критерии точности не превышающем Сумма Модель Модель Модель Модель Модель Модель Модель баллов №1 №2 №3 №4 №5 №6 №7 по всем моделям Рунге 2 3 2 3 3 2 2 Кутт Адамс 4 4 4 4 4 3 4 Башфорт Адамс 1 2 3 2 2 4 3 Моултон ФДН 3 1 1 1 1 1 1 Таблица 2.6 – Ранжирование методов при критерии точности не превышающем Сумма Модель Модель Модель Модель Модель Модель Модель баллов №1 №2 №3 №4 №5 №6 №7 по всем моделям Рунге Кутт 3 2 3 3 3 3 3 Адамс Башфорт 4 4 4 4 4 2 4 Адамс Моултон 1 3 2 2 2 4 2 ФДН 2 1 1 1 1 1 1 Ниже приведено сравнение методов решения систем ОДУ по точности при h=0,0001 и при tє[0, 10] (таблица 2.7, 2.8). Наиболее точным методом является метод Рунге-Кутта.

Таблица 2.7 – Сумма относительных ошибок Модель Модель Модель Модель №1 №2 №3 № Рунге-Кутт 0,00616 2,27E-08 0,02159 0, Адамс-Башфорт 0,01204 0,01967 0,04759 0, Адамс-Моултон 0,05095 0,04619 0,07516 3, ФДН 0,01162 0,01137 0,03644 1, Таблица 2.8 – Ранжирование методов по точности Модель Модель Модель Модель Сумма баллов №1 №2 №3 №4 по все моделям Рунге-Кутт 1 1 1 2 Адамс 3 3 3 1 Башфорт Адамс 4 4 4 4 Моултон ФДН 2 2 2 3 С увеличением размерности системы время, затраченное на ее решение, растет. В результате анализа графиков зависимости времени решения от количества уравнений можно сделать вывод, что это экспонента (Рисунки 2.3, 2.4). Система для сравнения методов по скорости выглядит следующим образом:

dP K i M 1 K p M 1 P1 K t P dt............................................

dPx K p M 1 ( Px 1 Px ) K t Px, 2 x dt dM K i M 1 K p M 1 Pk K t P dt k.............................................

dM x K t Px dt Рисунок 2.3 – Зависимость времени решения системы ОДУ от количества уравнений в системе ФДН время, с ФДН 100 200 300 400 Количество уравнений в системе Рисунок 2.4 – Зависимость времени решения системы ОДУ от количества уравнений в системе методом ФДН 2.4 Алгоритмы выбора длины шага Применение переменного шага интегрирования позволяет учитывать характер поведения решения и уменьшать общее число шагов, сохранив при этом требуемую точность приближенного решения. Таким способом могут быть снижены объем работы и машинное время и замедлен рост вычислительной погрешности [49].

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

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

Таким образом, критерий останова выглядит следующим образом (2.2), назовм его критерием останова по изменению функции:

y y предыдущее текущее J (2.2) 1 yтекущее Скорость изменения функции можно рассчитать по следующей формуле:

yi y y y i1 i1 i J (2.3) 2 yi y i На каждом шаге пересчитывать значения критерия для каждой функции требует значительных вычислительных затрат, поэтому обычно функцию принимают за константу не на одном шаге, а на отрезке.

Чаще всего применяется алгоритм выбора шага с помощью удвоения и деления шага пополам.

Пусть n+1 — оценка локальной погрешности метода на шаге h в точке (хn + h). Если pn 1 — наперед заданной погрешности, то считается, что не удовлетворяет заданной точности и нужно выбрать новый шаг, вдвое меньше прежнего h(1)=h/2. Вновь просчитывается значение в новой точке. Если новое значение погрешности, то точка и значение опять исключаются из рассмотрения, шаг снова делится пополам и вычисления повторяются. Так продолжается до тех пор, пока при какой-то величине шага оценка локальной погрешности n+1 не станет меньше или равной :

После этого считается, что решение дифференциального уравнения продолжено до точки Если же оценка локальной погрешности на шаге становится много меньше, то есть удовлетворяет неравенству |n+1| /M, где M — некоторая константа, то считается, что достигнута точность, превышающая заданную, и шаг интегрирования удваивается:

Если выполняется неравенство то считается, что полученное в точке xn+1 решение удовлетворяет заданной точности, и шаг интегрирования остается без изменения:. Дальнейшее решение уравнения производится из точки xn+1 с шагом hn+1=h(k).

Таким образом, обеспечивается выбор величины шага в зависимости от характера поведения решения дифференциального уравнения. Обычно полагается M = 2p, где p – порядок аппроксимации метода.

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

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

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

Описанные выше методы были апробированы на следующей билинейной динамической системе (2.4).

(2.4) Векторы выходных координат модели и матрицы коэффициентов приведены ниже [69]:

Для решения модели был разработан программный комплекс, написанный на объектно-ориентированном языке программирования С++.

В качестве эталона возьмм решение методом Рунге-Кутта четвртого порядка при шаге h = 0,00001. Ниже приведены эталонные решения модели (рисунки 2.5, 2.6).

Рисунок 2.5 – Решение модели при tє[0, 1] Рисунок 2.6 – Решение модели при tє[0, 10] Для расчта критерия, характеризующего точность решения модели, надо выполнить следующие шаги:

1. Берм по 100 контрольных точек с эталона и решения исследуемого метода.

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

.

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

, где n – порядок химической системы, k – количество контрольных точек.

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

Таблица 2.9 – Экспериментальные данные, полученные при использовании критерия останова, основанного на изменении функции Количество вычислений 480 000 2 080 000 4 640 Способ выбора шага Первый способ 53,0875% 58,6125% 0,6875% Второй способ 148,925% 0,8125% 0,15% Третий способ 9,4625% 0,5% 0,1125% Таблица 2.10 – Экспериментальные данные, полученные при использовании критерия останова, основанного на скорости изменения функции Количество вычислений 480 000 2 080 000 4 640 Способ выбора шага Первый способ 1,8125% 0,09% 0,025% Второй способ 3% 0,1425% 0,000763% Третий способ 0,04% 0,012% 0,0001% Так как временные затраты на вычисления также являются важной характеристикой при оценке различных методов, то произведем экспериментальные расчеты производительности при использовании различных критериев останова (таблица 2.11, 2.12).

Таблица 2.11 – Время, затраченное на вычисления, в мс при использовании критерия останова, основаного на изменении функции Количество вычислений 480 000 2 080 000 4 640 Способ выбора шага Первый способ 702 1623 Второй способ 702 1575 Третий способ 749 1670 Таблица 2.12 – Время, затраченное на вычисления, в мс при использовании критерия останова, основаного на скорости изменения функции Количество вычислений 480 000 2 080 000 4 640 Способ выбора шага Первый способ 718 1622 Второй способ 718 1747 Третий способ 718 1653 Как видно из экспериментальных данных, критерий, основанный на скорости изменения функции, показал немного более точное решение по сравнению с критерием останова, основанном на изменении функции. Однако обычно на практике используют комбинацию этих методов.

Первый способ выбора длины шага является достаточно грубым методом, поскольку длина шага может изменяться только в N-ое количество раз. При критерии останова, основаном на скорости изменения функции при количестве вычислений равном 480 000, получили не совсем гладкие кривые (рисунок 2.7):

Очевидным недостатком второго способа выбора динны шага является то, что функция может быть принятой за константу на слишком большом отрезке. В итоге при критерии останова, основаном на скорости изменения функции, при количестве вычислений равном 480 000, мы получили (рисунок 2.8).

Наилучшим методом является третий способ выбора длины шага. Он лишен выше описанных недостатков предыдущих способов. В итоге при критерии останова, основанного на скорости изменения функции при количестве вычислений равном 480 000, получили гладкие кривые (рисунок 2.9).

Рисунок 2.7 – Решение модели при использовании первого способа выбора длины шага Рисунок 2.8 – Решение модели при использовании второго способа выбора длины шага Рисунок 2.9 – Решение модели при использовании третьего способа выбора длины шага 2.5 Алгоритм переменного порядка и шага При решении задачи Коши для жестких систем ОДУ большой размерности в ряде случаев возникает необходимость применения алгоритмов на основе явных методов. Алгоритмы интегрирования на основе неявных или полуявных формул, как правило, используют обращение (декомпозицию с выбором главного элемента по строке или столбцу, а иногда и по всей матрице) матрицы Якоби, что в данном случае есть отдельная трудоемкая задача. Затраты на декомпозицию порядка N арифметических операций, где N есть размерность исходной системы. Кроме того, получение элементов матрицы Якоби и составление подпрограммы ее нахождения требуют от вычислителя больших затрат личного времени. В такой ситуации предпочтительнее применять алгоритмы на основе явных формул, если жесткость задачи позволяет за разумное время получить приближенное решение [86].

При одинаковом шаге явные методы в зависимости от порядка метода имеют разную устойчивость: чем выше порядок, тем ниже устойчивость метода (рисунок 2.10). Однако, чем выше порядок явного метода, тем выше точность вычислений. Логично предположить, что использование методов более низкого порядка на неустойчивых участках приведет к увеличению устойчивости метода по сравнению с методами более высоких порядков. А использование явных методов более высоких порядков на устойчивых участках приведет к повышению точности по сравнению с методами более низких порядков.

Описанный выше алгоритм был разработан программно с использованием методов Рунге-Кутта до 4-го порядка включительно и вошл в состав разработанного программного комплекса, описанного в главе 4.

Рисунок 2.10 – Область устойчивости явных методов Адамса Результаты вычислений представлены в таблице 4.1. Как видно из полученных данных, разработанный алгоритм позволил объединить достоинства явных методов Рунге-Кутта разных порядков. Полученный метод устойчив, как явный метод Рунге-Кутта первого порядка и при этом точен, как метод Рунге Кутта 4-го порядка.

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

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

Очевидно, что если потратить «освобожднное» время на дополнительные вычисления, то решение системы будет более точным.



Pages:   || 2 |
 





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

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