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

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

Новиков Антон Евгеньевич

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

01.01.07 - вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

2 2 ПАЯ 2014

Красноярск - 2014

005548681

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте вычислительного моделирования СО РАН (ИВМ СО РАН)

Научный руководитель: доктор физико-математических наук,

член-корреспондент РАН Шайдуров Владимир Викторович

Официальные оппоненты: Добронец Борис Станиславович

доктор физико-математических наук, профессор,

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

Задорин Александр Иванович

доктор физико-математических наук, профессор,

Омский филиал Федерального государственного бюджетного учреждения науки Института математики им. С.Л. Соболева Сибирского отделения Российской академии наук, заведующий лабораторией математического моделирования в механике

Ведущая организация: Федеральное государственное образовательное

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

Защита диссертации состоится 25 июня 2014 г. в 15 часов на заседании диссертационного совета Д 003.061.01 на базе Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики Сибирского отделения Российской академии наук (ИВМиМГ СО РАН) по адресу: г. Новосибирск, проспект академика Лаврентьева, 6.

С диссертацией можно ознакомиться в библиотеке ИВМиМГ СО РАН и на сайте www.sscc.ru.

Автореферат разослан 24 апреля 2014 г.

Ученый секретарь диссертационного совета

д.ф.-м.н. Рогазинский Сергей Валентинович

Общая характеристика работы

Актуальность темы обусловлена необходимостью эффективного численного решения жестких и нежестких систем обыкновенных дифференциальных уравнений большой размерности, возникающих при математическом моделировании динамики электрических сетей и электронных схем, биологических, химических и других процессов. Класс задач, описываемых жесткими системами, расширяется за счет учета большого числа факторов при построении математических моделей различных процессов. Это повышает требования к вычислительным алгоритмам. Современные методы решения жестких задач, как правило, на каждом шаге требуют обращение матрицы Якоби, что при большой размерности задачи определяет общие вычислительные затраты. С целью экономии во многих алгоритмах одна матрица Якоби применяется на нескольких шагах интегрирования. Эта задача решается достаточно просто в методах, в которых данная матрица используется в некотором итерационном процессе (многошаговые методы Адамса и Гира, неявные и полуявные методы типа Рунге-Кутты). В безытерационных схемах матрица Якоби включена в численную формулу, и поэтому возникают сложности с ее замораживанием. Однако такие методы обладают хорошими свойствами точности и устойчивости, просты с точки зрения реализации и, как следствие, привлекательны для вычислителей. Аналогом замораживания является применение в расчетах алгоритмов интегрирования на основе явных и I-устойчивых методов с автоматическим выбором численной формулы. Первые алгоритмы такого типа для многошаговых методов были построены Ре1го1с1 ЬЛ. (1982), для одношаговых численных схем - 8Ьатрте Ь.Б. и Новиков Е.А. (1984). В работах Ре1гоИ Ь.Я. и БИатрте Ь.Р. для выбора схемы применялась норма матрицы Якоби, которую нужно вычислять даже при расчетах по явной схеме. В работах Новикова Е.А. для выбора схемы применяется неравенство для контроля устойчивости без вычисления матрицы Якоби. По эффективности расчетов такие алгоритмы существенно превосходили существовавшие методы. Важно, что они распознавали, является ли задача жесткой или нет, и в зависимости от этого выбирали на каждом шаге эффективную численную формулу — явную или ¿-устойчивую. Поэтому построение новых эффективных алгоритмов переменного порядка и шага, а также переменной структуры, является актуальной задачей.

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

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

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

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

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

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

На защиту выносятся следующие результаты, соответствующие паспорту специальности 01.01.07 —вычислительная математика.

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

2. Явные методы первого порядка с расширенными областями устойчивости с неравенствами для контроля точности и устойчивости.

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

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

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

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

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

Теоретическая ценность. На основе методов Фельберга и Дорманда-Принса пятого, седьмого и восьмого порядков точности построены

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

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

Основные положения и результаты диссертации докладывались и обсуждались на Международных конференциях «Математические и информационные технологии» (Сербия, 2011), "Вычислительные и информационные технологии в науке, технике и образовании" (Павлодар, 2006), «Современные проблемы прикладной математики и механики: теория, эксперимент и практика» (Новосибирск, 2011), «Актуальные проблемы прикладной математики, информатики и механики» (Воронеж, 2009, 2010, 2011, 2012, 2013), «Студент и научно-технический прогресс» (Новосибирск, 2009, 2010, 2011, 2012), на Всероссийских конференциях «Имитационное моделирование. Теория и практика» (Санкт-Петербург, 2011), "Математическое моделирование развития северных территорий Российской Федерации" (Якутск, 2012), «Математическое моделирование и информационные технологии» (Красноярск, 2011, 2012, 2013); на семинарах Института вычислительного моделирования СО РАН и кафедры математического обеспечения дискретных устройств и систем Сибирского федерального университета. Работа поддержана грантами РФФИ (проекты 11-01-00106 и 14-01-00047).

Основные результаты диссертации опубликованы в 18 печатных работах, включая (в скобках в числителе указан общий объем этого типа публикации в печатных листах, в знаменателе - объем принадлежащий лично автору) 9 статей в периодических изданиях, рекомендованных ВАК (4.0/2.9) и 8 статей в трудах конференций (2.0/1.8).

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

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

СОДЕРЖАНИЕ ДИССЕРТАЦИИ

Во введении дан обзор работ по теме диссертации и приведено краткое описание ее содержания по главам.

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

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

В первом параграфе для численного решения задачи Коши

У = fit,у), Я'о) = Уо> 1о-{-1к (1)

изучается метод типа Рунге-Кутты вида

>Vl = Уп + РтхК + - + Рт А . (2)

где

ki=W(.t„,y„), ki=hf(t„ +^h,yn

3 3 9

К = ¥(к+-Ку„+—к, +—кг),

з 7v" g " 32 1 32 2

, ,,, 12 , 1932, 7200, 7296, ч

к, = hf(t +—h, v +-к.--к, +-к,), (3)

4 J^n 13 2197 1 219? 2 2197 37, v >

, ,/v , 439 ; о, 3680 , 845 , ч

К =hf(t+h,y„ +-к. -8к, +-к,--к.),

5 п, 216 1 2 513 3 4104 47

, ,/v 8 , 3544 / 1859 , 11 /ч K=hf(t+-h,y„--к, +2L--к, +-к,--кЛ.

в 2 27 2 2565 3 4104 4 40 5

При значениях коэффициентов

_ 16 _ 6656 _ 28561 ___£_ __2_

Рп ~ 135 ' Рп ' Аз " 12825' Ръ4 " 56430' Pss ~ 50' " 55 численная формула (2) совпадает с методом Фельберга и имеет пятый порядок точности. Известен другой набор коэффициентов

25 л 1408 2197 1

pn=-, р42= 0, Р«= —, Р»=Ш> л. = о,

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

где 11-11 — некоторая норма в Rw. В результате для контроля точности применяется неравенство г.п5 < е, где е - требуемая точность расчетов. Из результатов расчетов данным алгоритмом умеренно жестких задач следует, что на участке установления решения возникает большое количество возвратов из-за возникающей неустойчивости численной схемы.

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

ЧИСЛа Vn hÄnjпах

матрицы Якоби системы (1). Оценка получена на линейной задаче у —Ау, и вычисляется через стадии (3) по формуле

v„ = тах|(32&з -48к2 +I6kt)(/(к2 -кх).|/9.

Полученная оценка грубая, и поэтому она применяется как ограничитель на величину шага интегрирования. В результате шаг выбирается следующим образом. Учитывая, что е„ 5=0(А5) и v„=0(h), определяются два числа q и г по формулам q5e„ i=e и rv„= 3.6, где числом 3.6 ограничен интервал устойчивости метода Фельберга. Тогда прогнозируемый шаг /?„,, по точности и устойчивости выбирается по формуле

/|„+1 = тах[йя, mm(q,r)h„] . (4)

Данная формула позволяет стабилизировать размер шага на участке установления решения. Дополнительный контроль устойчивости приводит к повышению эффективности расчетов умеренно жестких задач примерно в полтора раза. Однако на участках установления точность вычислений значительно выше задаваемой точности. Это естественно, потому что старые ошибки подавляются за счет контроля устойчивости, а новые невелики за счет малости производных решения. В такой ситуации выгоднее считать методом более низкого порядка точности, но с более широкой областью устойчивости. На основе стадий (3) построена схема первого порядка точности с коэффициентами ри=0.41975960186956, р,2=0.44944365216575, р,3=0.1296419611922,jp14=0.1219923563523M0-2,^i5^-0.66250690732054-10"t, Pi6=0.11118997045939-10~5, интервал устойчивости которой расширен до 72 по действительной оси. Для контроля точности вычислений метода первого порядка применяется неравенство |1-2c2|x||£2-£i||<£, а ПРИ выборе величины шага интегрирования дополнительно проверяется |1-2c2|x||///[}vm)--ä:i||<£. Так как длина интервала устойчивости схемы первого порядка равна 72, то для ее контроля устойчивости применяется неравенство v„<72.

Методы первого и пятого порядков точности основаны на одних и тех же стадиях (3). Поэтому алгоритм переменного порядка и шага формулируется тривиально. При расчетах по методу пятого порядка нарушение неравенства v„<3.6 вызывает переход на схему первого порядка. При расчетах по методу первого порядка выполнение неравенства v„<3.6 вызывает переход на схему пятого порядка. При расчетах по методу первого порядка дополнительно контролируется устойчивость, а шаг выбирается по формуле вида (4). Так как интервал устойчивости схемы первого порядка примерно в 20 раз шире метода пятого порядка, то при расчетах алгоритмом переменного порядка и шага теоретическое повышение эффективности расчетов примерно в 20 раз по сравнению с исходным методом Фельберга.

Во втором и третьем параграфах аналогичные алгоритмы построены на основе стадий метода Фельберга седьмого порядка и метода Дорманда-Принса восьмого порядка точности. В четвертом параграфе приведены результаты расчетов. Через Fel5 и Fel7 обозначены алгоритмы интегрирования переменного шага на основе методов Фельберга пятого и седьмого порядков, а через DP8 — метод Дорманда-Принса восьмого порядка точности. Данные алгоритмы применительно к решению нежестких задач широко известны, их программные реализации входят практически во все мировые библиотеки. Здесь данные методы применяются для решения жестких задач. В качестве тестового примера выбрана простейшая математическая модель описания реакции Белоусова-Жаботинского (орегонатор). Задача является слишком жесткой для явных методов (коэффициент жесткости примерно и поэтому для ее решения

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

у[=П.П{уг-уху2 л-у, -8.375-10"V), У 2 = (у} ~ У2 ~ У ¡У2) / 77.27, = 0.1610>, - >>з), t е [0,300], h0 = 10~3, я(0) = 4, >-2(0) = 1.1, _у3(0) = 4. Через Fel5st, Fel7st и DP8st обозначены соответствующие алгоритмы интегрирования с дополнительным контролем устойчивости, а через Fel5vo, Fel7vo и DP8vo — алгоритмы интегрирования переменного порядка и шага. В качестве критерия эффективности выбрано число вычислений правой части задачи на интервале интегрирования. Время счета пропорционально данному критерию. Результаты расчетов приведены в табл. 1, а зависимость переменной yt от времени на рис. 1.

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

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

Таблица I

Вычислительные затраты

Точность расчетов lO"2 10"4 10"6 10~8

Fel5 15 694 434 15 691 105 15 704 622 15 692 418

Fel5st 12 855 329 12 871 206 12 890 292 12 934 168

Fel5vo 826 849 892 643 922 846 1 095 739

Fel7 38 429 365 38 429 235 38 436 138 38 461 306

Fel7st 19 836 063 19913 816 20 020 143 20 182 863

Fel7vo 778 253 830 494 1 046 225 1 574 532

DP8 27 999 975 27 995 836 27 993 163 28 008 372

DP8st 18 895 258 18 986 507 19 111 152 19 306 419

DP8vo 717 695 748 840 929 711 1 654 860

1ооооо

воооо

ROODO

4DOOO

20000

Рис. 1. Зависимость переменной у\ от времени (фрагмент)

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

Рассматривается задача Коши вида

/ = /00, У(1«) = Уо, Ч ■ (5)

Рассмотрение автономной задачи не снижает общности - введением дополнительной переменной неавтономную задачу можно привести к автономному виду, ¿-устойчивый (2,1)-метод второго порядка имеет вид

У„+1 = У» + аК + (! -<*)к2, D„kx = hf(y„), D„k2 =kt, (6)

где a=l-0.5V2, D„=E-ahA,„ E — единичная матрица, h — шаг интегрирования, An - матрица, представимая в виде An=f'n+hB„+0(hr), f',i=qf[y,)lду, В„ -некоторая матрица, не зависящая от шага интегрирования. Данное условие позволяет применять (6) с замораживанием как аналитической, так и численной матрицы Якоби. Неравенство для контроля точности имеет вид

Оценка максимального собственного числа м>„д=кЛ,и ,шх матрицы Якоби системы (5), необходимая для перехода на явную формулу, оценивается через ее норму по формуле w„$=h\\c}f{y)lcty\\.

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

У„+, = У„ +0.5(к, +кг), k,=hf{tn,yn), k2=hf(tn+h,yn+ki), (7) а неравенство для контроля точности - O^Hfe-AiH^f. Неравенство для контроля устойчивости схемы (7) построено с применением вспомогательной стадии кз=И/(уп+\). Оценка максимального собственного числа \vnj=hX,u гаах матрицы Якоби системы (5) вычисляется по формуле

wn,2 = 2max \(к} - к2)./(к2 -.

Интервал устойчивости схемы (7) приблизительно равен двум. Поэтому для ее контроля устойчивости применяется неравенство w„j2<2.

Метод первого порядка имеет вид

Уп+1=У„+( 7^+Л2)/8, (8)

где стадии к\ и fe описаны в (7). В качестве многочлена устойчивости метода (8) применен полином Чебышева степени 2. Поэтому область устойчивости (8) расширена до 8 по вещественной оси. Оценку максимального собственного числа w„_i=/?2„max матрицы Якоби системы (5) можно вычислить по формуле

= 8 ~ kl I kl ^ I'

Интервал устойчивости численной схемы (8) равен восьми. Поэтому для ее контроля устойчивости применяется неравенство w„,i<8.

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

в случае выполнения неравенства w„j<2. При расчетах по методу первого порядка наряду с точностью контролируется устойчивость w„ i<8.

В случае использования ¿-устойчивой схемы (6) формулировка алгоритма интегрирования также не вызывает трудностей. Нарушение неравенства "и>„д<8 вызывает переход на схему (6). Передача управления явным методам происходит в случае выполнения неравенства м>„>0<8, где оценка и^о вычислена через норму матрицы Якоб и системы (5). Численную формулу (6) без потери порядка точности можно применять с замораживанием матрицы Д,. В четвертом и пятом параграфах построены аналогичные алгоритмы на основе ¿-устойчивых и явных методов третьего и четвертого порядков.

В главе 4 приведены результаты сравнения эффективности построенных алгоритмов с методами Гира. Построенные алгоритмы переменной структуры обозначены соответственно ткгк2, шкгкЗ и ткгк4, а программа с реализацией метода Гира сПзоёе взята из библиотеки ЫейлЬ. Вычислительные затраты приведены в форме ¡ОД), где ¡Г - число вычислений правой части исходной задачи, у - количество декомпозиций матрицы Якоби. В первом параграфе описан алгоритм формирования дифференциальных уравнений химической кинетики. Во втором параграфе приведены результаты моделирования пиролиза этана, который в отсутствии кислорода описывается шестью стадиями при участии восьми реагентов. Соответствующая система состоит из 8 дифференциальных уравнений. В начале интервала интегрирования наблюдается переходный участок (сотые доли секунды), а затем происходит медленное установление. Коэффициент жесткости примерно 10'. Вычислительные затраты представлены в табл. 2, а поведение концентрации Н от времени приведено на рис. 2.

Таблица 2

Вычислительные затраты при моделировании пиролиза этана

Точность вычислений 1<Г2 ю-4 10"6

сИвойе 109(9) 129(9) 186(12)

ткгк2 61(4) 194(10) 3612(11)

ткгкЗ 36(3) 72(5) 179(12)

шкгк4 46(5) 63(6) 87(7)

Рис. 2. Зависимость концентрации Н от времени

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

Модель описывается пятью обратимыми и одной необратимой стадией при участии семи реагентов. Данная реакция протекает в изотермическом реакторе постоянного объема с обменом вещества, то есть соответствующая система дифференциальных уравнений состоит из семи уравнений. Коэффициент жесткости примерно 108. Вычислительные затраты представлены в табл. 3, а поведение концентрации [ВГ] от времени приведено на рис. 3.

Таблица 3

Вычислительные затраты при моделировании орегонатора

Точность вычислений ю-2 10Г* 10"6

ёЬоёе 7 806(542) 13 057(761) 27 205(1 524)

ткгк2 5 623(583) 13 749(989) 122 916(3 600)

ткгкЗ 7 566(595) 18 488(835) 103 665(7 130)

ткгк4 6 777(520) 13 157(629) 29 639(1 525)

1

н ц 1. Ц| к. х J ч, к

Рис. 3. Зависимость концентрации [ВЮ2] от времени

В четвертом параграфе приведены результаты моделирования проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма. Рассматривается система одномерных уравнений реакции-диффузии, которые возникают из химической реакции А+В—>С, где А - антитело с радиоактивной меткой, реагирующее с субстратом В - тканью, пораженной опухолью. При выводе уравнений предполагалось, что кинетика реакции описывается законом действующих масс, причем реагент А подвижен, тогда как реагент В неподвижен. После некоторых преобразований данная задача методом прямых сводится к задаче Коши для системы 400 обыкновенных дифференциальных уравнений. Коэффициент жесткости примерно 106. Вычислительные затраты представлены в табл. 4, а решение задачи проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма приведено на рис. 4.

В пятом параграфе приведены результаты моделирования кольцевого модулятора. Задача возникла из анализа электрических схем. Получая на входе низкочастотный сигнал и!п\ и высокочастотный сигнал иы2, кольцевой

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

Таблица 4

Затраты при решении задачи проникновения радиоактивных меток

Точность вычислений ю- 10"4 10 6

сИэоёе 25 358(62) 37 470(92) 64 021(157)

ткгк2 14 106(38) 53 749(189) 92 916(536)

ткгкЗ 23 749(59) 48 467(135) 73 665(177)

ткгк4 16 575(41) 31 157(87) 69 639(172)

Рис. 4. Решение задачи проникновения радиоактивных меток

Таблица 5

Вычислительные затраты при моделировании кольцевого модулятора

Точность вычислений 102 КГ* 10"6

с1Ьос1е 21 862(1 067) 28 867(1 267) 39 312(1 537)

ткгк2 30 135(1 782) 152 157(928) 184 973 (1 861)

ткгкЗ 21 971(734) 46 168(1 464) 48 609(1 895)

ткгк4 27 095(833) 28 604(1 117) 42 423(1 531)

Рис. 5. Зависимость выходного сигнала и2 от времени

Из анализа результатов расчетов всех четырех задач следует, что построенные алгоритмы переменного порядка, шага и переменной конфигурации при точности расчетов 1(Г2 и 10^ эффективнее программы dlsode. При точности вычислений Ю-6 на трех задачах эффективнее метод Гира. Это связано с тем, в dlsode реализованы методы более высокого порядка, чем в построенных алгоритмах.

Заключение содержит основные результаты и выводы диссертации.

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

Публикации в периодических изданиях, рекомендованных ВАК:

1. Новиков А.Е. Алгоритм переменного порядка и шага на основе стадий метода Дорманда-Принса восьмого порядка точности / А.Е. Новиков, Е.А. Новиков // Вычислительные методы и программирование. — 2007. — Т. 8. — С. 317-325.

2. Новиков А.Е. ¿-устойчивый (2,1)-метод решения жестких неавтономных задач / А.Е. Новиков, Е.А. Новиков // Вычислительные технологии. - 2008. - Т. 13. - Вестник КазНУ. - №3(58). - С. 477^82.

3. Новиков А.Е. Численное моделирование цикла цезия в верхней атмосфере ¿-устойчивым методом второго порядка точности / А.Е. Новиков, Е.А. Новиков // Вестник СибГАУ. - 2009. - Часть 1, №4(24). - С. 77-80.

4. Новиков А.Е. Численное решение жестких задач с небольшой точностью / А.Е. Новиков, Е.А. Новиков // Математическое моделирование. — 2010. — Т. 22, №1,-С. 46-56.

5. Novikov А.Е. Numerical Integration of Stiff Systems with Low Accuracy / A.E. Novikov, E.A. Novikov // Mathematical Models and Computer Simulations. — 2010. - Vol. 2, No. 4. - P. 443—452.

6. Новиков A.E. Численное моделирование пиролиза этана (2,1)-методом решения жестких неавтономных задач / Л.В. Кнауб, А.Е. Новиков, Ю.А. Шитов // Вестник КрасГАУ. - 2010. - №1. - С. 22-27.

7. Кнауб Л.В. Семейство (т,1)-методов решения жестких линейных задач / Л.В. Кнауб, А.Е. Новиков // Вестник ИжГТУ. - 2010. - №4(48). - С. 152155.

8. Новиков А.Е. Комбинированный алгоритм третьего порядка для решения жестких задач / А.Е. Новиков, Е.А. Новиков // Вычислительные технологии. - 2011. - Т. 16, №6. - С. 54-68.

9. Новиков А.Е. Алгоритм интегрирования переменной конфигурации на основе явно-неявных схем четвертого порядка / А.Е. Новиков, В.В. Шайдуров // Вестник Бурятского государственного университета. — 2012. -Выпуск9.-С. 111-116.

10. Новиков А.Е. Контроль устойчивости метода Ческино второго порядка точности / А.Е. Новиков, Е.А. Новиков, Л.В. Кнауб // Вестник Бурятского государственного университета.-2013.-Выпуск 9.-С. 111-116.

Публикации в трудах конференций:

11. Новиков А.Е. Замораживание матрицы Якоби в (2,2)-методе решения жестких систем // Международная конференция «Актуальные проблемы прикладной математики, информатики и механики». - 2009. - Воронеж, 22-24 июня. - С. 87-89.

12. Новиков А.Е. Модифицированный метод Дорманда-Принса переменного порядка и шага / А.Е. Новиков // Международная конференция «Актуальные проблемы прикладной математики, информатики и механики»,- 2010.— Воронеж, 20-22 сентября.- С. 274-277.

13. Новиков А.Е. Неоднородный метод третьего порядка для решения жестких задач / А.Е. Новиков // ХЫХ Международная научная студенческая конференция "Студент и научно-технический прогресс". -2011,- Новосибирск, 16-20 апреля. - С. 218-219.

14. Новиков А.Е. Согласование областей устойчивости явных методов / А.Е. Новиков, Е.А. Новиков // Международная конференция «Математические и информационные технологии». - 2011. - Сербия -Черногория, 27 августа - 5 сентября. http://conf.nsc.ru/MIT-2011/геротче\у/47120.

15. Новиков А.Е. Алгоритм интегрирования переменной структуры для решения жестких задач / А.Е. Новиков, Е.А. Новиков // Пятая Всероссийская конференция по имитационному моделированию и его применению в науке и промышленности «Имитационное моделирование. Теория и практика» (ИММОД-2011). - 2011. - Санкт-Петербург, 19-21 октября.-С. 195-200.

16. Новиков А.Е. Численный метод третьего порядка для решения жестких задач / А.Е. Новиков // Сборник трудов международной конференции «Актуальные проблемы прикладной математики, информатики и механики». - 2011. - Воронеж, 26-28 сентября. - С. 280 -282.

17. Новиков А.Е. Применение явных трехстадийных методов для моделирования кинетики химических реакций / А.Е. Новиков // III Всероссийская научная конференция студентов, аспирантов, молодых ученых и специалистов "Математическое моделирование развития северных территорий Российской Федерации". - 2012. - Якутск, 21 - 28 мая. С. 127-130.

18. Новиков А.Е. Явно-неявный алгоритм четвертого порядка для решения жестких задач / А.Е. Новиков // Международная конференция «Актуальные проблемы прикладной математики, информатики и механики». - 2012. - Воронеж, 17-19 сентября. - С. 274—278.

Подписано в печать 24.04.2014 г. Формат 60x84/16. Объем 1 п.л. Тираж 100 экз. Отпечатано на ризографе ИВМ СО РАН 660036, г. Красноярск, Академгородок, 50/44

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ УЧРЕЖДЕНИЕ НАУКИ ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОГО МОДЕЛИРОВАНИЯ СИБИРСКОГО ОТДЕЛЕНИЯ РОССИЙСКОЙ АКАДЕМИИ НАУК

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

04201455707

Новиков Антон Евгеньевич

РАЗРАБОТКА АЛГОРИТМОВ ПЕРЕМЕННОЙ СТРУКТУРЫ ДЛЯ

РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ

01.01.07 - вычислительная математика

ДИССЕРТАЦИЯ

на соискание ученой степени кандидата физико - математических наук

Научный руководитель:

доктор физико-математических наук, член-корреспондент РАН Шайдуров Владимир Викторович

Красноярск - 2014

Оглавление

Введение..................................................................................................................3

Глава 1. Контроль точности и устойчивости одношаговых методов.......16

1.1. Основные определения..................................................................................16

1.2. Контроль точности вычислений....................................................................21

1.3. Контроль устойчивости..................................................................................24

1.4. Реализация явных методов с контролем точности и устойчивости..........27

Глава 2. Алгоритмы интегрирования переменного порядка и шага.......29

2.1. Алгоритм на основе схемы Фельберга пятого порядка..............................29

2.2. Алгоритм на основе схемы Фельберга седьмого порядка..........................34

2.3. Алгоритм на основе схемы Дорманда-Принса восьмого порядка............43

2.4. Результаты расчетов.......................................................................................52

Глава 3. Комбинированные алгоритмы интегрирования..........................56

3.1. Методы типа Розенброка...............................................................................56

3.2. Класс (ш,к)-методов решения жестких задач..............................................60

3.3. Численное решение жестких задач с небольшой точностью.....................62

3.3.1. L-устойчивый (2,1)-метод.......................................................................62

3.3.2. Явный метод второго порядка................................................................65

3.3.3. Явный метод первого порядка................................................................68

3.3.4. Алгоритм интегрирования переменной структуры..............................69

3.4. Комбинированный алгоритм третьего порядка...........................................71

3.4.1 L-устойчивый (3,2)-метод.........................................................................71

3.4.2. Явный метод третьего порядка...............................................................76

3.4.3. Явный метод первого порядка................................................................82

3.4.4. Алгоритм интегрирования переменной структуры..............................83

3.4.5. Замораживание матрицы Якоби в (3,2)-методе....................................84

3.5. Комбинированный алгоритм четвертого порядка.......................................87

3.5.1. Максимальный порядок точности (ш,2)-методов................................88

3.5.2. L-устойчивый (4,2)-метод........................................................................90

3.5.3. Метод Мерсона с контролем точности и устойчивости......................94

3.5.4. Явный метод первого порядка................................................................95

3.4.5. Алгоритм интегрирования переменной структуры..............................98

Глава 4. Результаты расчетов практических задач...................................100

4.1. Дифференциальные уравнения химической кинетики.............................101

4.2. Пиролиз этана................................................................................................103

4.3. Модифицированный орегонатор.................................................................104

4.4. Проникновение помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма.............................................106

4.5. Кольцевой модулятор...................................................................................108

Заключение.........................................................................................................112

Список литературы..........................................................................................113

Введение

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

У = Rt,у), y(t0) = y0,t0<t<tk, (1)

где у и /- вещественные TV-мерные вектор-функции, t - скалярная величина. Учет большого числа факторов при построении математических моделей приводит к расширению класса задач, описываемых жесткими системами [3, 5— 6, 9-25]. Основные тенденции при построении численных методов связаны с решением систем большой размерности [1-3, 6-7, 12-14, 17—21]. Сложность практических задач приводит к возрастающим требованиям к вычислительным алгоритмам, поэтому проблема создания эффективных численных методов решения задачи Коши для жестких систем большой размерности является актуальной задачей. При построении эффективных алгоритмов интегрирования требуется разрешить ряд вопросов [1-5, 18-21]. Нужно выбрать методы, соответствующие решаемым задачам. Здесь рассматриваются одношаговые безытерационные схемы вида [5]

Уп+1=Уп + к<РА*п,Уп^), п = 0,1,2,.... (2)

На некоторых задачах они имеют преимущество перед многошаговыми схемами [3-5, И, 17-18, 20, 23], которые приводят к осреднению решения (срезание экстремумов), что при моделировании некоторых динамических объектов делает их неприемлемыми [4]. Если правая часть задачи (1) разрывная, то применение многошаговых методов малоэффективно [4]. Требование безытерационности (2) позволяет оценить затраты на шаг до проведения расчетов и упрощает программную реализацию [1-2, 23].

В последнее время появилось огромное количество методов

интегрирования жестких задач [1-5, 13-14, 16, 19-31]. Для перехода от идеи метода к его алгоритмической реализации нужно решить важный круг вопросов, связанных с изменением величины шага интегрирования и оценкой точности получаемых численных результатов [1-5, 13-14, 19-21, 32-47]. Современные способы управления шагом основаны, как правило, на контроле точности численной схемы [1-5]. Такой подход представляется наиболее естественным, поскольку основным критерием при проведении практических расчетов является точность нахождения решения. Многие алгоритмы интегрирования для выбора величины шага интегрирования используют оценку локальной ошибки или погрешности аппроксимации. Это оправдано тем, что если на каждом шаге контролировать некоторый минимальный уровень локальной ошибки, то глобальная ошибка будет ограничена. Можно выделить три практических способа оценки данной ошибки [1-2].

Классическим является способ, основанный на экстраполяционной формуле Ричардсона. Его еще называют правилом Рунге, и он заключается в следующем [5, 48-49]. В каждой "сеточной точке интервала интегрирования решение вычисляется с шагом к и 0.5/г, а искомая оценка определяется через разность приближений к решению. Недостатком способа является необходимость дважды вычислять решение в каждой точке, что приводит к значительному увеличению вычислительных затрат.

Более дешевым является многошаговый способ [1-2, 50]. Он заключается в том, что одношаговой формуле р-то порядка точности в соответствие ставится многошаговая схема (р+1)-го порядка. Затем данная схема преобразуется таким образом, чтобы после подстановки в нее приближений (2) получилась оценка локальной ошибки одношагового метода. Недостатком данного способа является многошаговость оценки, что приводит ко всем недостаткам многошаговых методов.

В последнее время все большую популярность получает способ оценки локальной ошибки с помощью вложенных методов [5, 51]. Приближение к решению в каждой точке вычисляется двумя методами р-то и (р+1)-го порядков

точности вида (2). Затем локальная ошибка метода р-то порядка оценивается через разность полученных приближений. Обычно такой способ используется тогда, когда расчеты по методу р-то порядка не требуют дополнительных вычислений правой части и матрицы Якоби исходной задачи. По вычислительным затратам на шаг такой способ лежит между оценкой ошибки с помощью правила Рунге и многошаговой оценкой. Однако, по отношению к многошаговой оценке, в ней при вычислении ошибки используется информация только с данного шага, что повышает ее надежность. Данный способ хорошо зарекомендовал себя при решении многих задач [1-5, 19-21, 52-58] и ниже будет использоваться здесь.

Использование оценки локальной ошибки для контроля точности вычислений в ряде случаев приводит к успеху. Однако с целью повышения надежности необходимо найти оценку глобальной ошибки [1—5, 41-47]. Наиболее известный способ определения данной ошибки основан на предположении о линейном характере накопления глобальной ошибки из локальных ошибок [5, 59-60]. В результате для контроля точности предлагается использовать неравенство \\дп\\<£к/^-Ь), где 3„ — оценка локальной ошибки, ||-|| - некоторая норма в Яы, е - требуемая точность расчетов. Очевидно, что предположение о линейном характере накопления является достаточно грубым. Иногда используется другое неравенство ||<5Л||<фСр|1/р [61-63], которое получено в предположении, что при интегрировании устойчивыми методами вклад начальных возмущений убывает по геометрической прогрессии.

Интенсивное исследование неявных методов началось после работы Дальквиста [64], в которой было введено понятие ^-устойчивости. Это понятие привело к рассмотрению неявных методов типа Рунге-Кутты [65-70]

Известна теорема - для каждого т существует неявная т-стадийная схема (3) порядка 2т. Аналогичная теорема для явных методов типа Рунге-Кутты отсутствует. Методы максимального порядка могут быть только А-

\

]=

У

(3)

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

Несмотря на хорошие свойства точности и устойчивости, практическое использование неявных методов типа Рунге-Кутты является ограниченным вследствие больших вычислительных затрат на шаг интегрирования. При реализации (3) требуется применение итерационного процесса. Метод простой итерации не эффективен при решении жестких задач, так как он приводит к такому же ограничению на размер шага, что и явный метод. Поэтому необходимо использование метода Ньютона-Рафсона или каких-либо его модификаций. Это приводит к необходимости обращения матрицы размерности {тИхтЩ, где т - число стадий в (3), N - размерность системы (1). Сокращение вычислительных затрат достигается за счет использования одной матрицы на нескольких шагах интегрирования. Замораживание матрицы становится возможным благодаря тому, что это не влияет на порядок точности численной схемы, а определяет только скорость сходимости итераций. Поэтому необходимость в ее пересчете возникает при значительном замедлении сходимости итерационного процесса.

Трудности с реализацией неявных методов типа Рунге-Кутты привели к поискам более простых их модификаций. Был рассмотрен класс полуявных формул типа Рунге-Кутты, для которых в (3) имеет место Ду=0 при /</. В этом случае итерационная матрица является блочно-диагональной, причем число блоков совпадает с числом стадий, а размерность каждого блока с размерностью вектора решения. В результате, вместо обращения матрицы размерности {тЫхтЫ), надо обратить т матриц размерности N каждая. Дальнейшего сокращения вычислительных затрат можно добиться, если положить равными все /?„-, 1 <1<т, и аппроксимировать все диагональные матрицы одной. Тогда на шаге требуется обратить только одну матрицу размерности (тУхТУ). В этом случае порядок точности (т+2) не может быть достигнут ни для какого т-стадийного полуявного метода при /?ц =.. .= ртт.

Розенброком были предложены два метода второго и третьего порядков

точности [31]. В дальнейшем такие схемы стали называть методами типа Розенброка, и они имеют вид

где Е - единичная матрица, д/!ду - матрица Якоби задачи (1), рь аг, у у, 1 <1<т, 1</</-1, - числовые коэффициенты. Для упрощения данные методы записаны для автономной системы. Наиболее эффективные реализации методов типа Розенброка возникают тогда, когда все аь 1 <1<т, равны между собой, а уу=0. В этом случае на каждом временном шаге требуется вычисление и обращение всего лишь одной матрицы размерности (тУхЛО. Для улучшения свойств точности таких методов Ваннером предложено вычислять стадии по следующим формулам [71]

Данная модификация получила название ЯС^-методов. Методы типа Розенброка и ЯОАУ-методы относятся к одношаговым безытерационным методам. Они существенно более просты с точки зрения реализации, чем алгоритмы на основе численных формул, в которых стадии вычисляются с использованием итераций. Однако в данных методах матрица Якоби влияет на порядок точности численной схемы. Поэтому их недостатки в основном связаны с проблемой замораживания матрицы £>„, декомпозиция которой при больших N фактически определяет общие вычислительные затраты. Если вопрос об использовании одной матрицы на нескольких шагах оставить нерешенным, то в случае использования безытерационных методов нужно ограничиваться задачами небольшой размерности. На основе схемы типа Розенброка второго порядка точности построен алгоритм интегрирования с замораживанием матрицы Якоби [72]. Однако при построении такого алгоритма на основе формулы третьего порядка возникают принципиальные трудности. Применение 110\У-методов проблемы не решает.

Еще одним важным требованием к современным алгоритмам

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

Для явных методов шаг интегрирования к ограничен неравенством ^И-тах!^ [3-5], где Лтю1 есть максимальное собственное число матрицы Якоби системы (1), а положительная постоянная £> связана с размером области устойчивости. Так как для многих жестких задач длина интервала интегрирования значительно превышает величину £>/|Л.тах|, то интегрирование при условии к\Хта^<Л) оказывается непосильным для современных ЭВМ. В последнее время в связи с построением явных методов с расширенными областями устойчивости их возможности трактуются более широко [73—89].

Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и Ь-устойчивых методов с автоматическим выбором численной схемы [89—93]. В этом случае эффективность алгоритма может быть повышена за счет расчета переходного участка, соответствующего максимальному собственному числу матрицы Якоби, явным методом. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости. Применение таких комбинированных алгоритмов полностью не снимает проблему замораживания матрицы Якоби, потому что явным методом

можно просчитать, вообще говоря, только погранслойное решение, соответствующее максимальному собственному числу [3—4].

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

Актуальность темы обусловлена необходимостью эффективного численного решения жестких и нежестких систем обыкновенных дифференциальных уравнений большой размерности. Современные методы решения жестких задач, как правило, на каждом шаге требуют обращение матрицы Якоби, что при большой размерности задачи определяет общие вычислительные затраты. С целью экономии во многих алгоритмах одна матрица Якоби применяется на нескольких шагах интегрирования. Аналогом замораживания является применение в расчетах алгоритмов интегрирования на основе явных и 1-устойчивых методов с автоматическим выбором численной формулы. Первые алгоритмы такого типа для многошаговых методов были построены Ре1го1с1 Ь.Я. (1982), для одношаговых численных схем — БЬашрте ЬЛ\ и Новиков Е.А. (1984). В работах Ре1го1с1 Ь.Я. и ЗЬашрте ЬЛ7. для выбора схемы применялась норма матрицы Якоби, которую нужно вычислять даже при расчетах по явной схеме. В работах Новикова Е.А. для выбора схемы применяется неравенство для контроля устойчивости без вычисления матрицы Якоби. По эффективности расчетов такие алгоритмы существенно превосходили существовавшие методы. Важно, что они распознавали, является ли задача жесткой или нет, и в зависимости от этого выбирали на каждом шаге эффективную численную формулу - �