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

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

V

Г

17

Кнауб Людмила Владимировна

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

05.13.18 - математическое моделирование, численные методы и комплексы

программ

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

1 4 ^ДЯ 2СС9

Красноярск 2009

003469217

Работа выполнена в Институте вычислительного моделирования СО РАН, г.Красноярск и Федеральном государственном образовательном учреждении высшего профессионального образования «Сибирский федеральный университет», г.Красноярск.

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

профессор

Новиков Евгений Александрович

Официальные оппоненты: доктор физико-математических наук,

профессор

Добронец Борис Станиславович

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

Черных Геннадий Георгиевич

Ведущая организация: Тюменский государственный университет

Защита диссертации состоится «29» мая 2009 г. в 14 часов в ауд. УЖ-115 на заседании диссертационного совета ДМ 212.099.06 ФГОУ ВПО «Сибирский федеральный университет» по адресу: 660074, г. Красноярск, ул. Киренского, д. 26 б.

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

Ваш отзыв заверенный печатью предприятия просим направлять по адресу: 660074, г. Красноярск, ул. Киренского, 26 б, учёному секретарю диссертационного совета ДМ 212.099.06. тел./факс. (3912)912-235.

Автореферат разослан "27" апреля 2009 .

Учёный секретарь

диссертационного совета Р.Ю.Царев

о

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

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

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

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

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

• Разработаны новые методы с расширенными областями устойчивости.

• Для данных методов построены неравенства для контроля точности

вычислений и устойчивости численной схемы.

• Созданы алгоритмы переменного порядка и шага.

• Разработан комплекс программ для решения задач средней жесткости.

Теоретическая значимость

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

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

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

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

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

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

Результаты диссертации докладывались:

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

Семинар кафедры защиты информации Тюменского государственного университета.

Семинары отдела вычислительной математики Института вычислительного моделирования СО РАН (Красноярск).

Семинары кафедры математического обеспечения дискретных устройств и систем (МОДУС) Сибирского федерального университета.

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

На защиту выносится:

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

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

• Комплекс программ с программными реализациями построенных методов интегрирования.

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

Общая характеристика диссертации. Работа состоит из введения, пяти глав, заключения, библиографического списка из 105 наименований и приложения. Общий объем работы составляет 122 стр., приложения на 1стр.

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

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

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

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

y' = f(t,y\ yQ0) = y0, t0<,t<tk, (1)

у' = КУ), У(0 = Уо, (2)

применяются явные методы типа Рунге-Кутта следующего вида, соответственно

т /-1

Упл = Уп + Ти РтЛ. К = ¥(1п + а,Куп + £ ДЛ,), (3)

т (-1

р-А. =+Е АЛА (4)

М /=1

где ртП а„ Ру,\<1<т,\<]<1-\, - коэффициенты, определяющие свойства точности и устойчивости (3) и (4), И - шаг интегрирования, - стадии метода, Для упрощения выкладок все исследования проводились для методов (4). Однако все построенные численные схемы можно использовать для решения неавтономных систем. Для этого достаточно в (3) подставить а„1<г<т, определенные следующим образом,

о,= 0, = 2<1<т.

№

Во втором параграфе построен двухстадийный метод

Уп+1 =Уп + + к2), к, = = + £,)> (5)

неравенство для контроля точности которого имеет вид

0.5|| ||<£-, (6)

где || -1| - некоторая норма в й", £ - требуемая точность вычислений. В третьем параграфе разработан метод третьего порядка

2,

К = ¥(У„), К = ¥(У„ + ^А), (7)

неравенство для контроля точности которого имеет вид

Ьк1-2к2 + к3\\<Е. (8)

о

В четвертом параграфе для пятистадийного метода Мерсона

1 , 2, 1 ,

(?)

к5=к/{Уп+1-кх-иг+2к<):

построено неравенство для контроля точности

2^-9*3 + ^4-^ 5 А

(10)

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

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

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

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

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

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

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

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

В первом параграфе получена оценка максимального собственного числа уп 2 = ИЯп тгх матрицы Якоби для метода (5). Неравенство для контроля

устойчивости имеет вид

уП2 = 2тах(|к3-к'2 \1\к!1-к[\)<2, (11)

где числом 2 ограничен интервал устойчивости схемы (5), а кг = И/(уп^) есть вспомогательная стадия, необходимая для оценки максимального собственного числа. Вычисление к3 к дополнительным вычислительным затратам не приводит, потому что къ совпадает со стадией к{ для следующего шага.

Во втором параграфе получена оценка максимального собственного числа у„ з = ИЯп шах матрицы Якоби для метода (7). Неравенство для контроля устойчивости имеет вид

у„ з = 0.5тах(| к[ - 2к[ + к'3 | /1 к'2 - к[ |) < 2.5. (12)

В третьем параграфе получена оценка максимального собственного числа \п_4 = тах матрицы Якоби для метода (9). Неравенство для контроля устойчивости имеет вид

у„ 4 = 6шах(| к\ - к[ | /1 к'2 - к\ |< 3.5. (13)

Для всех методов сформулирован алгоритм интегрирования с контролем точности вычислений и устойчивости численной формулы. Прогнозируемый шаг Ип+1 вычисляется по следующей формуле

^тах^.ттС/Л/,1')], (14)

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

точности, К' - прогнозируемый шаг по устойчивости. Из формулы (14) следует, что контроль устойчивости применяется как ограничитель на рост шага, что позволяет сгладить последствия грубости оценок собственных чисел матрицы Якоби. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг В результате для выбора шага и предлагается формула (14). Данная формула позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Собственно говоря, именно наличие данного участка существенно ограничивает возможности применения явных методов для решения жестких задач.

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

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

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

В первом параграфе построен алгоритм интегрирования на основе численных формул первого и второго порядков точности. Метод первого порядка имеет вид

= Уп К = ¥(У„), К = ¥{У„ + Ю- (15)

Для контроля точности метода (15) применяется неравенство 3|| к2 -к{ || /8< е. Область устойчивости метода первого порядка расширена

до 8 по вещественной оси. Если шаг ограничен требованием устойчивости, то переход с (5) на (15) может приводить к повышению эффективности расчетов в 4 раза. Переключение со схемы (5) на (15) происходит в случае нарушения неравенства у„ 2 < 2. При расчетах по методу первого порядка контролируется неравенство у„, < 8, где числом 8 ограничен интервал устойчивости метода (15), , = 8шах(| к'г-к[\1\к[- к[ |). Шаг

интегрирования выбирается по формуле (14). Обратное переключение со схемы (15) на (5) происходит в случае выполнения неравенства у„д < 2.

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

517, 208, 4 ,

У, 4.1 = У„ +-к. +-л, +-я,,

729 1 729 2 729 3

=¥(У„~к1+2к2),

Для контроля точности метода (16) применяется неравенство 19||А2-^||/27<£. Область устойчивости метода первого порядка расширена до 18 по вещественной оси. Если шаг ограничен требованием устойчивости, то переход с (7) на (16) может приводить к повышению эффективности расчетов в примерно в 7 раз. Переключение со схемы (7) на (16) происходит в случае нарушения неравенства у„ 3 <2.5. При расчетах по методу первого порядка контролируется неравенство уп3 <18, где числом 18 ограничен интервал устойчивости метода (16). Шаг интегрирования выбирается по формуле (14). Обратное переключение со схемы (16) на (7) происходит в случае выполнения неравенства ул3 < 2.5.

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

+ (17)

где стадии &(,1<г'<т, определены формулами (9), р, = 5,248365568• 10"1, рг =3.260928-10"', 1.395154944-10"1, Л =9.5158272-10'3,

р5 =3.93216 -10"5. Для контроля точности метода (17) применяется неравенство НА^-^ргч Область устойчивости метода первого порядка расширена до 50 по вещественной оси. Если шаг ограничен требованием устойчивости, то переход с (9) на (17) может приводить к повышению эффективности расчетов в примерно в 14 раз. Переключение со схемы (9) на

(17) происходит в случае нарушения неравенства vr:A <3.5. При расчетах по методу первого порядка контролируется неравенство уп4 <50, где числом 50 ограничен интервал устойчивости метода (17). Шаг интегрирования выбирается по формуле (14). Обратное переключение со схемы (17) на (9) происходит в случае выполнения неравенства у„ 4 < 50.

В четвертом параграфе приведены результаты численного эксперимента. Вычислительные затраты приведены в табл. 1-3. В табл. 1-3 используются обозначения:

ЯК2 - метод второго порядка с контролем точности вычислений; ИС28Т - метод второго порядка с контролем точности и устойчивости; ИК2РР - алгоритм переменного порядка и шага;

ИКЗ - метод третьего порядка с контролем точности вычислений; МСЗБТ - метод третьего порядка с контролем точности и устойчивости; ШСЗРР - алгоритм переменного порядка и шага;

ЯК4 - метод четвертого порядка с контролем точности вычислений; ИК48Т - метод четвертого порядка с контролем точности и устойчивости; ЯК4РР - алгоритм переменного порядка и шага;

Заметим, что программа ЯК4 совпадает с известным методом Мерсона, программная реализация которого МЕК^ОИ входит во все известные библиотеки мира. Из анализа табл. 1 - 3 можно сделать вывод о том, что применение алгоритмов переменного порядка и шага для решения задач умеренной жесткости приводит к существенному повышению эффективности и, в частности, эффективнее программы МЕЯБОМ примерно в 15 раз (см. строки 7 и 9 табл. 1-3).

Таблица 1

Суммарные затраты с точностью s = 10~2

Метод Число шагов Число Число вычислений

возвратов функции/

RK2 3 922 747 1 272 100

RK2ST 4 101 759 1 894 8 205 814

RK2PP 1 104 327 2 726 2 211 387

RK3 3 049 717 807 768 10 764 687

RK3ST 3 115 679 8 638 9 364 313

RK3PP 458 596 1864 1 377 692

RK4 2 160 775 809 885 14 051 457

RK4ST 2 361 823 13 197 12 180 374

RK4PP 170948 742 855 935

Таблица 2

Суммарные затраты с точностью е = КГ*

Метод Число шагов Число Число вычислений

возвратов функции/

RK2 3 925 209 1 276 748 9 127 209

RK2ST 4 104 091 3 075 8 211 336

RK2PP 1 053 099 3 753 2 109 958

RK3 3 050 308 807 164 10 765 252

RK3ST 3 115 834 11 344 9 370 190

RK3PP 466 671 2 687 1 402 782

RK4 2 161 031 814 393 14 062 734

RK4ST 2 387 296 16 059 12 159 808

RK4PP 171 302 892 857 536

Таблица 3

Суммарные затраты с точностью е = КГ6

Метод Число шагов Число Число вычислений

возвратов функции/

RK2 3 947 998 1 277 186 9 173 189

RK2ST 4 128 568 6 312 8 263 473

RK2PP 1 148 263 5 723 2 302 256

RK3 3 052 669 807 618 10 773 243

RK3ST 3 117 083 18812 9 388 873

RK.3PP 483 414 9 165 1 459 692

RK4 3 010 085 809 729 14 045 063

RK4ST 2 480 079 22 151 12 489 006

RK4PP 182 630 2 891 917 576

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

На рис. 1. используются следующие обозначения:

ИК_ООЕ - управляющая программа;

ИК2РР - программная реализация алгоритма интегрирования переменного порядка и шага на основе явной двухстадийной схемы типа Рунге-Кутта;

RK2ST1 - программная реализация алгоритма интегрирования переменного шага на основе явной двухстадийной схемы типа Рунге-Кутта первого порядка точности;

RK2ST - программная реализация алгоритма интегрирования переменного шага на основе явной двухстадийной схемы типа Рунге-Кутта второго порядка точности;

RK2ST1

RK3ST1

RK ODE

RK3ST

RK4ST1

RP(n, t, у, f)

Рис. 1. Структурная схема комплекса программ ЯК ОБЕ

ШСЗРР - программная реализация алгоритма интегрирования переменного порядка и шага на основе явной трехстадийной схемы типа Рунге-Кутта;

ШС38Т1 - программная реализация алгоритма интегрирования переменного шага на основе явной трехстадийной схемы типа Рунге-Кутта первого порядка точности;

ЯКЗБТ - программная реализация алгоритма интегрирования переменного шага на основе явной трехстадийной схемы типа Рунге-Кутта третьего порядка точности;

МС4РР - программная реализация алгоритма интегрирования переменного порядка и шага на основе явной пятистадийной схемы типа Рунге-Кутта;

ШС48Т1 - программная реализация алгоритма интегрирования переменного шага на основе явной трехстадийной схемы типа Рунге-Кутта первого порядка точности;

ЮС45Т - программная реализация алгоритма интегрирования переменного шага на основе явной пятистадийной схемы типа Рунге-Кутта четвертого порядка точности;

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

Отметим, что программы ЯК2РР, 11КЗРР и ЯК4РР сами по себе являются комплексами программ, в каждый из которых включено по пять алгоритмов переменного шага (расчеты методом первого порядка без контроля устойчивости; методом первого порядка с контролем устойчивости; методом высокого порядка без контроля устойчивости; методом высокого порядка с контролем устойчивости; расчеты алгоритмом переменного порядка и шага).

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

А + У^±Х + Р, ¿, = 0.084, Х + Г&2 Р, к2 = 4-Ю8, А + Х^± 2Ж, *з=2-10\ С + Жт=>Х + г, =1.3 105, 2 Х^±А + Р, к5= 4-Ю7, г->С + 0.462Г, к6 = 0.65.

кА = Ю\ А:_2 = 5-10"5, к_3 = 2 • 107, к_4 = 2.4 • 107, к_5 = 4-10~",

где кь -5<1<6, - константы скоростей прямых (с положительными индексами) и обратных (с отрицательными индексами) стадий. В данной реакции участвуют 7 частиц, имеющие следующие обозначения:

А=ВгСГз, С=М(п), Р=НОВг, УГ= Вг02, Х=НВЮ2, У~Вг~, 2=М{п+\).

В этих обозначениях М(п) - ион металла катализатора, М(п+\) - окисленная форма этого иона. Обозначим концентрации реагентов следующим образом: с\=[ВгСГз], с2=[ВГ[, с3= [М(п)], с4= [НВЮг], с5=[НОВг], с6= [Вг02], с7=[Л/(и+1)]. Данная реакция протекает в изотермическом реакторе постоянного объема с обменом вещества, то есть ей соответствует система из семи дифференциальных уравнений вида

С = АУг+^(Ср-С),

где С=(сь с2, ..., с/ - вектор концентраций реагентов, А -стехиометрическая матрица, К=(У|, ..., - вектор скоростей стадий,

Ср={ср 1, ........ср1)т - вектор концентраций реагентов на входе в реактор, & -

время пребывания смеси в реакторе, 0~г/и, г - объем реактора, и - объемная скорость течения смеси через реактор.

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

са = -V, -у2 + 0.462у6 + (ср2 -с2)1в, С'3=-У4+У6+(СР3-С3)/9,

с4 = П "уз + ~2у5 +(СР4 С'5=У1 + 2У1+У!+(Ср5-С5)/0, С^2г3-У4+(Ср6-С6)/0, С7=Ъ-У6+(ср7-с7)/в,

где скорости V], у2.....стадий определяются по формулам

У1 = КС£г ~ ^-1С4С5 ' У2 = ^2С2С4 - ^-2С5 > У3 = ^3С1С4 ~ >

У4 = ^4С3С6 ~ ^-4С4С7 > У5 ^5С4 ~ ^-5С1С5 > У6 = КС1 '

Концентрации реагентов на входе в реактор принимают значения сР1=0.14, сР2=0.15М0"5, сру= 0.125-10-3, ^=^5=^=^7=0, причем 0=125.5. Начальные значения концентраций реагентов следующие: С1=0.1387, £Г2=0-1534-10"6, с3=0.1176-10"3, с4=0.3165-10-7, с5=0.1956-10"3, с6=0.5814-10'6, с7=0.631-10"5. На данной задаче КК._СЮЕ эффективнее МЕИБСЖ примерно в 11 раз. Зависимость концентрации [Вг~] от времени приведена на рис. 2.

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

ди _ дги _ ^ 5у

= -^-кт, — = -кт, (18)

которые возникают из химической реакции А+В—>С с константой скорости реакции к, где А- антитело с радиоактивной меткой, реагирующее с субстратом Л - тканью, пораженной опухолью. Концентрации А и В обозначены через и и V соответственно. При выводе уравнений (18) предполагалось, что кинетика реакции описывается законом действующих масс, причем реагент А подвижен, тогда как реагент В неподвижен.

Изучается полубесконечная пластина, внутри которой субстрат В равномерно распределен. Реагент А, попадая на поверхность пластины, начинает проникать в нее. Для моделирования проникновения уравнения (18) рассматриваются в полосе 8т~{{х,1)'. 0<х<оо, О<КТ} с начальными и(х, 0)=0, \{х, 0)=у0, х>0, и граничными и(0, ?)= Ф(0< 0<КТ , условиями, где у0 константа. После некоторых преобразований данная задача методом прямых сводится к задаче Коши для системы обыкновенных дифференциальных уравнений

~ = Яиу),у{0) = г, уеЯ2\0<1<20, (19)

м

где Ы- задаваемый пользователем параметр. Функция / определяется формулами

Г -а Угм-Уъ-з , п Уц->-2Угн+Угм у

^ 1 2А£ ^ (лО2

где

«у = 2(Ж -1 Ус2, р, = ЦЦ -1 Ус2,1<;<М,А£ = \/М, = Ш У2дч1 = Угы-V 8 е К™> 8 = (О,У0,О,У0, ...,О,У0)г.

Функция ф(?) = 2 при 0</<5 и ф(Г) = 0 при 5 </<20, то есть ф{г) имеет разрыв первого рода в точке ?=5. В качестве параметров к, у0 и с использовались ¿=100, у0=1 и с=4. Расчеты проводились при 200, то есть система (19) состоит из 400 уравнений. Задача о нахождении разрыва функции ф(*) при /=5 возлагалась на алгоритм управления шагом. Решение задачи (19) приведено на рис. 3.

Интегрирование задачи (19) с точностью £ = 10"2 проводилось на промежутке [0,20] с начальным шагом 10"5. Суммарные вычислительные затраты приведены в табл. 4. Расчеты проводились с применением

управляющей программы ЮСОВЕ. Для наглядности в табл. 4 приведены идентификаторы входящих в ИКОСЕ методов.

х-axis

Рис. 3. Решение задачи (19)

Таблица 4

Суммарные вычислительные затраты

Метод Число Число Число вычислений

шагов возвратов функции f

RK2 89 629 29 069 208 328 _

RK2ST 88 868 951 178 688

RK2PP 24 659 34 49 353

RK3 69 806 18 205 245 828

RK3ST 58 381 2 299 179 741

RK3PP 11 894 43 35 726

RK4 49 312 18 118 319 033

RK4ST 48 532 1 076 246 965

RK4PP 4 951 19 24 791

Напомним, что алгоритм ШС4 совпадает с программой МЕЯЗОИ, в которой реализован метод Мерсона. Из анализа таблицы 4 можно сделать вывод о том, что разработанный комплекс программ эффективнее известного метода Мерсона почти в 13 раз. Фактическая точность расчетов не хуже задаваемой точности.

В приложении приведены тестовые примеры.

Основные результаты.

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

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

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

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

4. Создан комплекс программ ЮСОВЕ, в состав которого включены программные реализации разработанных алгоритмов интегрирования.

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

Публикации по теме диссертации в журналах из списка ВАК:

1. Кнауб JI.B., Лаевский Ю.М., Новиков Е.А. Алгоритм интегрирования переменного порядка и шага на основе явного двухстадийного метода Рунге-Кутты // Новосибирск: СибЖВМ, т.10, №2,2007. - с. 177 - 185.

2. Кнауб Л.В., Новиков Е.А. Алгоритм интегрирования на основе явного трехстадийного метода Рунге-Кутта // Вестник КрасГАУ, №3,2009. - с. 49-54.

3. Кнауб Л.В., Новиков Е.А. Численное моделирование орегонатора двухстадийным методом типа Рунге-Кутта // Вестник КрасГАУ, №4, 2009.-с. 16-21.

4. Кнауб Л.В., Новиков Е.А. Применение явного трехстадийного метода типа Рунге-Кутта для численного моделирования задач химической кинетики // Вестник СибГАУ, №1(22), часть 1,2009. - с. 77-80.

5. Кнауб Л.В., Новиков Е.А. Алгоритм интегрирования с контролем точности и устойчивости явного трехстадийного метода типа Рунге-Кутта // Системы управления и информационные технологии, №1(35), 2009.-с. 20-24.

Публикации в других изданиях:

6. Кнауб Л.В., Новиков Е.А. Контроль устойчивости явного двухстадийного метода типа Рунге-Кутта // Труды международной конференции "Вычислительные и информационные технологии в науке, технике и образовании", т.2, Павлодар, 2006. - с. 71-76.

7. Кнауб Л.В., Новиков Е.А. Алгоритм интегрирования на основе явного двухстадийного метода типа Рунге-Кутта // Вестник КрасГАУ: Ресурсосберегающие технологии в с/х, №4, 2007. - с. 146-150.

Кнауб Людмила Владимировна Разработка эффективных явных методов и комплекса программ для решения задач химической кинетики умеренной жесткости Автореф. дисс. на соискание ученой степени кандидата физико-математических наук Подписано в печать 20.04.2009. Заказ ЗЪЗ Формат 60x90/16. Усл. печ. л. 1,2. Тираж 100 экз. Отпечатано в ИПК СФУ 660074, Красноярск, ул. Киренского, 28

Введение.

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

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

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

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

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

Глава 2. Алгоритмы с контролем точности вычислений.

2.1. Методы типа Рунге-Кутта.

2.2. Схемы второго порядка точности.

2.3. Схемы третьего порядка точности.

2.4. Алгоритм интегрирования на основе метода Мерсона.

2.5. Анализ результатов расчетов.

Глава 3. Алгоритмы с контролем устойчивости численной схемы.

3.1. Схема второго порядка точности.

3.2. Схема третьего порядка точности.

3.3. Схема четвертого порядка точности.

3.4. Анализ результатов расчетов.

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

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

4.2. Алгоритм интегрирования на основе трехстадийной схемы.

4.3. Алгоритм с применением стадий Рунге - Кутта - Мерсона.

4.4. Анализ результатов расчетов.

Глава 5. Комплекс программ RKODE и результаты моделирования практических задач.

5.1. Комплекс программ RKODE для решения жестких задач.

5.2. Численное моделирование реакции Белоусова-Жаботинского, дающей сложный предельный цикл.

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

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

У = fit,У), = h^t^h- С1)

Класс задач, описываемых жесткими системами, расширяется, так как учитывается все большее число факторов при построении математических моделей физических процессов. Поэтому возникает необходимость решения жестких задач все более высокой размерности. Это, в свою очередь, повышает требования к вычислительным алгоритмам. Современные методы решения жестких задач, как правило, на каждом шаге требуют обращение матрицы Якоби, что при достаточно большой размерности задачи определяет общие вычислительные затраты. Известные явные методы, в которых матрица Якоби не применяется, в основном не приспособлены для решения задач даже умеренной жесткости по двум причинам. Во-первых, области устойчивости явных методов малы, что приводит к обременительным ограничениям на величину шага интегрирования. Во-вторых, на участке установления решения шаг раскачивается из-за противоречивости требований точности вычислений и устойчивости численной схемы. Поэтому построение новых эффективных явных методов с расширенными областями устойчивости и контролем устойчивости численных схем, а также алгоритмов переменного порядка и шага, является актушНшг®Ьсщддаей.выбрать методы, соответствующие классу решаемых задач. Здесь в основу алгоритмов интегрирования положены явные схемы типа Рунге-Кутта, которые можно записать в виде

Уп+1=Уп + к(РА*п,Уп,Ь)> (2) где п — текущая точка интегрирования, h - шаг интегрирования, (pf — заданная вектор-функция, зависящая от правой части / задачи (1). Они обладают определенными преимуществами перед многошаговыми формулами. В частности, многошаговые методы приводят к осреднению решения ("срезание экстремумов"), что при моделировании некоторых динамических объектов делает их неприемлемыми. Обзор работ, посвященных численному решению задачи (1) многошаговыми методами, содержится в [19, 41, 44, 46, 56].

Существует большое количество методов интегрирования жестких систем. Важным! является круг вопросов, связанных с изменением шага интегрирования и оценкой точности получаемых численных результатов, что и делает метод экономичным и надежным. При проведении практических расчетов основным критерием является точность нахождения решения. Поэтому способы управления шагом основаны, как правило, на контроле точности численной схемы. Многие алгоритмы интегрирования при выборе величины шага используют оценку локальной ошибки или, что то же самое, погрешности аппроксимации, потому что если на каждом шаге контролировать некоторый минимальный уровень локальной ошибки, то глобальная ошибка будет ограничена. В настоящее время можно выделить три практических способа оценки данной ошибки ([41], с. 59-65).

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

V05l'-vh v05h-vh

-^- + 0(hp+2), -+ 0(hp+2)

1 — 1~p 2P - I где p - порядок точности метода, у''7 и у^51' - приближения к решению с шагом h и 0.5/г, S* и S^5ph - соответствующие локальные ошибки. Но этот способ приводит к значительному увеличению вычислительных затрат, т.к. необходимо дважды вычислять решение в каждой точке.

Более дешевым является многошаговый способ [57]. Он заключается в том, что одношаговой формуле р—го порядка точности в соответствие 4 ставится многошаговая схема (р +1) —го порядка

Е [мы - щт^))]=o(kp+2). i=0

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

8Я.Р = (Z Д )~'Е - ЩЯУпЛ (3)

0 /=0

Недостатком этого способа является многошаговость оценки.

В последнее время наиболее популярной является оценка локальной ошибки с помощью вложенных методов. Приближение к решению в каждой точке вычисляется двумя методами вида (2) р—то и (р +1)-го порядков точности, а затем локальная ошибка метода р -го порядка оценивается через разность полученных приближений дпр = упр+х ~Упр- Такой способ используется, когда для вычислений по методу р-го порядка не требуется дополнительных вычислений правой части и матрицы Якоби задачи (1). Следует отметить оперативность и относительную дешевизну оценки локальной ошибки с помощью вложенных методов. По затратам на шаг она лежит между оценкой ошибки с помощью экстраполяции Ричардсона и многошаговой оценкой. В то же время, по отношению к многошаговой оценке, в ней при вычислении ошибки используется информация только с данного шага, что повышает ее надежность. Данный способ успешно применялся в [41, 61, 63, 76, 85-86] и ниже будет использоваться здесь.

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

SJ<£h/(tk-t0), (4) где 8п — оценка локальной ошибки, ||-|| - некоторая норма в RN, б — требуемая точность расчетов. Такой способ успешно применялся в [56]. Очевидно, что предположение о линейном характере накопления является достаточно грубым. В [3] используется другое неравенство

Up

SJ\<£ s!cp

5) где р — порядок точности метода, а ср — некоторая вычисляемая постоянная. Это неравенство получено в предположении, что при интегрировании устойчивыми методами вклад начальных возмущений убывает по мере продвижения по сетке. Обоснование (5) для линейной скалярной задачи (1) приведено в ([35], с.43-45). Еще один способ оценки основан на интегрировании дополнительной линейной системы дифференциальных уравнений, описывающей поведение главного члена глобальной ошибки (см., например, [44], с. 40-45). Однако это связано с вычислением матрицы Якоби, которая в явных методах типа Рунге-Кутта не используется, и дополнительными затратами на интегрирование. Поэтому такой способ применяется достаточно редко.

При решении ряда задач L -устойчивыми методами возникает проблема с обращением матрицы Якоби. Поэтому при численном исследовании некоторых жестких задач все большее внимание привлекают явные методы (см. библ. [36]), которые не нуждаются в вычислении данной матрицы. Если жесткость задачи не слишком велика, то они будут предпочтительнее. Появление многопроцессорных вычислительных систем позволяет взглянуть иначе на явные методы, которые легко распараллеливаются [21]. Две основные причины, которые приводят к трудностям при использовании явных методов для решения жестких задач: а) противоречие между точностью и устойчивостью численной схемы на участке установления. Следствием этого является раскачивание шага интегрирования, что в ряде случаев заканчивается аварийной остановкой вычислений. Этого недостатка можно избежать, например, предложенным в [24] способом контроля устойчивости.

Ь) области устойчивости известных численных схем слишком малы. Работы [11, 13, 17, 21-22, 24-29, 31-37, 47-55, 59-60, 64-67, 69-72, 73-80, 87-94] посвящены вопросам построения явных методов с расширенными областями устойчивости.

Расширение области устойчивости связано с ростом вычислительных затрат на шаг интегрирования. Поэтому, если шаг ограничен по точности, такие схемы будут малоэффективны. Однако если шаг ограничен в силу устойчивости, что имеет место на участке установления, то за счет применения численных схем с расширенными областями устойчивости можно значительно повысить эффективность алгоритма интегрирования [17, 21, 28-29, 33, 36, 74, 90-92, 94]. В качестве критерия выбора численной формулы можно использовать неравенство для контроля устойчивости [24].

Очевидно, что за счет контроля устойчивости и использования численных схем с расширенными областями устойчивости можно только расширить границы применимости явных формул. Для данных методов шаг h ограничен в силу неравенства h | Лтах |< D, где Лтах есть максимальное собственное число матрицы Якоби системы (1), а положительная постоянная D связана с размером области устойчивости. Так как для многих жестких задач длина интервала интегрирования значительно превышает величину DI | Ятах |, то условие h | Хтах |< D является весьма и весьма обремКшитселмиБОДов, основанный на дробно — рациональном представлении приближенного решения исследуется в [73], а в [40] рассмотрены методы, основанные на использовании аппарата цепных и ветвящихся дробей. Еще один подход к построению вычислительных алгоритмов заключается в конструировании численных схем, учитывающих специфику исходной задачи. Здесь можно выделить схемы экспоненциальной подгонки [16], и методы, основанные на обращении главной части дифференциального оператора [8-9]. В [15] предложен новый численный метод интегрирования жестких систем, в основе которого лежит принцип последовательной фильтрации членов, соответствующих наибольшим собственным значениям матрицы Якоби системы (1). Он позволяет без потери устойчивости увеличить шаг интегрирования даже в случае простейших численных схем. В [2] рассмотрены вопросы реализации методов интегрирования на Фортране.

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

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

Основные результаты опубликованы в журналах по списку ВАК [99103].

Заключение

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

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

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

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

4. Создан комплекс программ RKODE, в состав которого включены программные реализации разработанных алгоритмов интегрирования.

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

1. Арушанян О.Б., Залеткин С.Ф., Захаров А.Ю., Калиткин Н.Н. О тестировании программ решения обыкновенных дифференциальных уравнений. / Препринт 139, М.: ИПМ АН СССР, 1983.

2. Арушанян О.Б., Залеткин С.Ф. Численное решение обыкновенных дифференциальных уравнений на Фортране. М.: Издательство Московского университета, 1990. 1975, с. 214-220.

3. Артемьев С.С., Демидов Г.В. Алгоритм переменного порядка и шага для численного решения жестких систем обыкновенных дифференциальных уравнений. //ДАН СССР, 1978, т. 238, 3, с. 517 520.

4. Артемьев С.С., Демидов Г.В. Алгоритм переменного порядка и шага для численного решения жестких систем обыкновенных дифференциальных уравнений. / Препринт 45, Новосибирск: ВЦ СО АН СССР, 1978.

5. Бахвалов Н.С. Численные методы. М.: Наука, 1975, т. 1,

6. Бобков В.В. Об одном семействе нелинейных разностных схем. // Дифференциальные уравнения, 1977, 11, с. 2086-2078.

7. Бобков В.В. Новые явные А устойчивые методы численного решения дифференциальных уравнений. // Дифференциальные уравнения, 1978, 12, с. 2249 -2252.

8. Бобков В.В.Об одном способе построения методов численного решения дифференциальных уравнений. // Дифференциальные уравнения, 1983, 7, с. 1115 1122.

9. Бобков В.В., Мандрик П.А., Репников В.И. Об одном классе разностных схем для дифференциальных уравнений. // Вестник Белорусского ун. та, 1985, Сер.1, Физ. матем. и мех. 3, с. 3 1 - 34.

10. Виноград Р.Э. Об одном критерии неустойчивости в смысле Ляпунова решений линейной системы обыкновенных дифференциальных уравнений. // ДАН СССР, 1952, т. 84, с. 201 204.

11. Деккер К., Вервер Я. Устойчивость методов Рунге — Кутта для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988.

12. Демидов Г.В., Новиков Е.А. Экономичный алгоритм интегрирования нежестких систем обыкновенных дифференциальных уравнений. // Числ. методы мат. физики, Новосибирск: ВЦ СО АН СССР, 1979, с. 69-83.

13. Демидов Г.В., Новиков Е.А. О контроле точности явных формул типа Рунге — Кутта второго и третьего порядков аппроксимации с помощью формул более низкого порядка. // Числ. методы мех. сплошной среды, Новосибирск, 1984, т. 15, 6, с. 59 -74.

14. Демидов Г.В., Новиков Е.А. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений. // Числ. мет. мех. сплошной среды, Новосибирск, 1985, т. 16, 1, с. 27-39.

15. Демирчян К.С., Ракитский Ю.В. О фильтрации составляющих с большими производными в динамических системах. / Препринт 3, М.: ИВТАН СССР, 1984.

16. Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методырешения задач с пограничным слоем. М.: Мир, 1983.

17. Дуракова В.К., Новиков В.А., Новиков Е.А. Явные методы типа Рунге Кутта первого порядка точности с заданным размером интервала устойчивости. // ЖВМ и МФ, 1988, т. 28, 4, с. 603 - 607.

18. Захаров А.Ю. Некоторые результаты сравнения эффективности методов решения систем обыкновенных дифференциальных уравнений. / Препринт 125, М.: ИПМ АН СССР, 1979.

19. Захаров А.Ю., Кальянова Н.А., Капуста В.О., Шульмина Т.П. О программах, комплексах и пакетах программ для решения обыкновенных дифференциальных уравнений. / Препринт 160, М.: ИПМ АН СССР, 1984.

20. Кальянова Н, А., Захаров А.Ю., Маркачев Ю.Е. LSODA пакет программ для численного решения жестких и нежестких систем обыкновенных дифференциальных уравнений. / Инструкция, М.: ИПМ АН СССР, 1988.

21. Лебедев В.И. Явные разностные схемы с переменными шагами по времени для решения жестких систем уравнений. / Препринт 177, М.: ОВМ АН СССР, 1987.

22. Локуциевский В.О., Локуциевский О.В. Применение чебышевских параметров для численного решения некоторых эволюционных задач. / Препринт 98, М.: ИПМ АН СССР, 1984.

23. Новиков Е.А. Некоторые эффективные алгоритмы численного решения задачи Коши для систем обыкновенных дифференциальных уравнений. / Дисс. канд. физ. мат. наук, Новосибирск: ВЦ СО АН СССР, 1983.

24. Новиков В.А., Новиков Е.А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных дифференциальных уравнений.//ДАН СССР, 1984, т. 277. 5, с. 1058- 1062.

25. Новиков Е.А. Построение алгоритма интегрирования жестких систем дифференциальных уравнений на неоднородных схемах. // ДАН СССР, 1984, т. 278, 2, с. 272-275.

26. Новиков В.А., Новиков Е.А. Об алгоритме переменной структуры на основе явных формул типа Рунге Кутта первого и второго порядков точности. / Препринт 112, Новосибирск: ВЦ СО АН СССР, 1985. '

27. Новиков В.А., Новиков Е.А. Явные методы для решения жестких систем обыкновенных дифференциальных уравнений. / Препринт 629, Новосибирск: ВЦ СО АН СССР, 1985.

28. Новиков Е.А., Шитов Ю.А. Исследование (т,к) — методов решения жестких систем с одним и двумя вычислениями правой части. / Препринт 15, Красноярск: ВЦ СО АН СССР, 1987.

29. Новиков В.А., Новиков Е.А. О построении явных методов типа Рунге

30. Кутта с расширенными областями устойчивости./ Препринт 9, Красноярск: ВЦ СО АН СССР, 1988.

31. Новиков В.А., Новиков Е.А. Численное конструирование областей устойчивости явных методов. / Препринт 15, Красноярск: ВЦ СО АН СССР, 1988

32. Новиков Е.А., Дуракова В.К. Алгоритм переменного порядка и шага на основе явной шестистадийной схемы типа Рунге Кутта. / Препринт 11, Красноярск: ВЦ СО АН СССР, 1988.

33. Новиков В.А., Новиков Е.А., Шокин Ю.И. Численное решение обыкновенных дифференциальных уравнений с небольшой точностью. // Вопросы качественной теории диф. уравнений, Новосибирск: Наука, 1988. с. 29-35.

34. Новиков В.А., Новиков Е.А. Одношаговые методы интегрирования задачи Коши для обыкновенных дифференциальных уравнений. / Красноярск: КГУ, 1989, ч.1, 1989.

35. Новиков В.А., Новиков Е.А., Дуракова В.К. Автоматизация построения на ЭВМ явных методов первого и второго порядков точности с адаптивной областью устойчивости. / Препринт 17 и 18. Красноярск: ВЦ СО АН СССР, 1989.

36. Новиков Е.А. Одношаговые безытерационные методы решения жестких систем. / Дисс. доктора физ. мат. наук, Новосибирск, 1991. — Р&]. Петровский И:Г. Лекции по теории обыкновенных дифференциальных уравнений. М.: Наука, 1970.

37. Ракитский Е.В., Устинов С.М., Черноруцкий И.Г. Численные методы решения жестких систем. М.: Наука, 1979.

38. Слоневский Р.В. Новые дробно рациональные численные методы решения жестких систем. // Числ. решение ОДЕ, М.: ИПМ АН СССР, 1988, с. 124-138

39. Современные численные методы решения обыкновенных дифференциальных уравнений. / Под ред. Дж. Холла и Дж. Уатта, М.: Мир, 1979.42. • Шилов Г.Е.Математический анализ (функции нескольких вещественных переменных). М.: Наука, 1972.

40. Широбоков Н.В. К определению жестких дифференциальных задач. //ЖВМиМФ, 1984, т. 24, 4, с. 599-601.

41. Штеттер X. Анализ методов дискретизации для обыкновенных дифференциальных уравнений. М.: Мир, 1978.

42. Федоренко Р.П. О регулярных системах обыкновенных дифференциальных уравнений. // ДАН СССР, 1983, т. 273, 6, с. 1318 -143692. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990.

43. Ansorge R., Tornig W. Zur stabilitat des Nystromschen verfahren. // Z. Angew. Math. Mech., 1960, 40, p. 568 570.

44. Bakker M. Analytic aspects of a minimax problem (Dutch). / Report TN 62, Amsterdam: Mathematisch Centrum, 1971.

45. Beentjes P.A., Dekker K.A. 5th order, 6th stage Runge Kutta formula with optimal stability boundary(Dutch). / Report NR 27. Amsterdam: Mathematisch Centrum, 1972.

46. Beentjes P.A. Some special formulas of the England class of fifth order Runge — Kutta schemes./ Report NW14/74, Amsterdam: Mathematisch Centrum, 1974.

47. Bettis D.G. Numerical integration of products of Fourier and ordinary polynomials.//Numer. Math., 1970, 14, p. 421 -434.

48. Bettis D.G. Runge Kutta algorithms for oscillatory problems. // Z. Angew. Math. Phys., 1979, 30, p. 699 - 704.

49. Burrage K. Order and stability properties of explicit multivalue methods. //Appl. Numer. Math., 1985, 1, p. 363 -379.

50. Butcher J.C. On the convergence of numerical solutions to ordinary differential equations. //Math. Сотр., 1966, 20, p. 1 10.

51. Dalquist G. A special stability problem for linear multistep methods. // BIT, 1963, v. 3, p. 23 43.

52. Dekker K. Generalized Runge Kutta methods for coupled systems of hyperbolic differential equations. // J. Сотр. Appl. Math., §6^|7,Bs]pke^K.-S2t&lbility of linear multistep methods on the imaginary axis. //BIT, 1981, 21, p. 66-69.

53. England R. Error estimates for Runge Kutta type solutions to systems of ODE's. // Comput. J., 1969, 12, p. 166 - 169.

54. Gear C.W. Hybrid methods for initial value problems in ordinary differential equations. // SIAM J. Numer. Anal., 1965, 2, p. 69 86.

55. Gear C.W. The automatic integration of stiff ordinary differential equations. //Infor. Proc., 1969, p. 187 193.

56. Gentzsch W., Schluter A. Ueber ein Einschrittverfahren mit zyklischer Schrittweitenanderung zur lozung parabolischer Differentialgleichungen. //ZAMM, 1978, 58, p. t415 -t416.

57. Henrici P. Discrete variable methods in ordinary differential equations. New York London: John Wiley & Sons., Inc. 1962.

58. Kinnmark I.P.E., Gray W.G. Fourth order accurate one - step integration methods with large imaginary stability limits. // Numer. Math, for Partial Differential Equations, 1986, 2. p. 63 - 70.

59. Lambert J. D. Computational methods in ordinary differential equations, Wiley, New York, 1973.

60. Lapidus L., Seinfeld J. H. Numerical solution of ordinary differential equation, Academic Press, New York, 1971.

61. Novikov V.A., Novikov E.A. On the accuracy and stability control of one step methods of integration of ordinary differential equations. //In.: Proc. BAIL-mConf., Bool Press, 1984, p. 81-93.

62. Novikov V.A., Novikov E.A. Explicit methods for stiff systems of ordinary differential equations. //In.: Proc. П Conf.: Numerical analysis and applications held in Benin City, Nigeria, Bool Press, 1986.

63. Novikov V.A., Novikov E.A. Explicit methods of Runge Kutta type with adaptive stability Region. //In.: Proc. BAIL - V Conf., Bool Press, 1988, p. 269-276.

64. Novikov E.A. Application of explicit Runge Kutta methods to solve stiff ODE's. // Advances in Modeling & Analysis, A, AMSE Press, v. 16, 1, 1992, p. 23 -35.

65. Richardson L.F. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stress in a massory dam. // Philos. Trans. Roy. Soc., London, 1910, ser. A, 210, p. 307 357.

66. Shintani H. On a one step method of order 4. // J. Sei. Hiroshima Univ. 1966. Ser. A - 1 Math., 30, p. 91 - 107.

67. Shintani H. Two step processes by oneOstep methods of order 3 and of order 4. // J. Sei. Hiroshima Univ., 1966, Ser. A-l Math., 30, p. PB] - l^&n. der Houwen P.J. One - step methods for linear initialvalueproblems. //ZAMM, 1971, 51, p. t59 -t60.

68. Van der Houwen P.J. Explicit Runge Kutta methods with increased stability boundaries. //Numer. Math., 1972, 20, p. 149 - 164.

69. Van der Houwen P.J. Construction of integration formulas for initial value problems. North - Holl and Amsterdam, 1977.

70. Van der Houwen P.J. Stabilized Runge Kutta methods for second - order differential equations without first derivatives. // SIAM J. Numer. Anal., 1979, 16, p. 523 - 537.

71. Van der Houwen P.J. Modified Nystrom methods for semi — discrete hyperbolic differential equations. // SIAM J. Numer. Anal. 1981, 18, p. 1081 1097.

72. Van der Houwen P.J., Sommeijer B.P. Predictor corrector methods with improved absolute stability regions. // IMA J. Numer. J^|l.ElrftS§iiS,¥-^7Hi4B7r.E. Comparing numerical methods for the solutions of systems of ODE's //BIT, №15, 1975 . - p. 10 - 48.

73. Корзухин M. Д., Жаботинский A. M. Математическое моделирование химических и экологических автоколебательных систем. — М.: Наука, 1965

74. Showalter К., Noyes R.M., Bar-Eli К. A Modified Oregonator Model Exhibiting Com plicated Limit Cycle Behavior in a Flow System // J. Chem. Phys., 69, 1978. - p. 2514 - 2524

75. Mazzia F., Magherini C. Test Set for Initial Value Problem Solvers //Department of Mathematics, University of Bari and INdAM, Research Unit of Bari, release 2.4, 2008.

76. Кнауб Л.В., Лаевский Ю.М., Новиков E.A. Алгоритм интегрирования переменного порядка и шага на основе явного двухстадийного метода Рунге-Кутты // Новосибирск: СибЖВМ, т. 10, №2, 2007.-е. 177 185. *

77. Кнауб Л.В., Новиков Е.А. Алгоритм интегрирования на основе явного трехстадийного метода Рунге-Кутта // Вестник КрасГАУ, №3, 2009. с. 4954.

78. Кнауб Л.В., Новиков Е.А. Численное моделирование орегонатора двухстадийным методом типа Рунге-Кутта // Вестник КрасГАУ, №4, 2009. -с. 16-21.

79. Кнауб Л.В., Новиков Е.А. Применение явного трехстадийного метода типа Рунге-Кутта для численного моделирования задач химической кинетики // Вестник СибГАУ, №1(22), часть 1, 2009. с. 77-80.

80. Кнауб Л.В., Новиков Е.А. Алгоритм интегрирования с контролем точности и устойчивости явного трехстадийного метода типа Рунге-Кутта // Системы управления и информационные технологии, №1(35), 2009. с. 2024.

81. Кнауб Л.В., Новиков Е.А. Контроль устойчивости явного двухстадийного метода типа Рунге-Кутта // Труды международной конф. "Вычислительные и информационные технологии в науке, технике и образовании", т.2, Павлодар, 2006. с. 71-76.

82. Кнауб Л.В., Новиков Е.А. Алгоритм интегрирования на основе явного двухстадийного метода типа Рунге-Кутта // Вестник КрасГАУ: Ресурсосберегающие технологии в с/х, №4, 2007. с. 146-150.