Актуальность проблемы. Отличительными чертами живых биологических систем, несмотря на их сложность, являются высокая упорядоченность и эффективная компактная сборка. В биологии и медицине часто встречаются примеры колебаний, которые возникают при самых разных обстоятельствах, и периоды которых варьируются от нескольких секунд до часов, дней и даже недель [1, 2]. Период этих колебаний может быть связан с периодическими изменениями условий жизни на Земле — смена времен года, смена дня и ночи. Существуют и другие геофизические ритмы — солнечные, лунные, связанные с периодами атмосферных явлений [3]. Но многие периодические процессы имеют частоту изменения, не связанную очевидным образом с внешними геокосмическими циклами [4]. Это так называемые «биологические часы» различной природы, начиная от колебаний биомакромолекул, биохимических колебаний, вплоть до популяционных волн [5]. Внутриклеточные колебания задают эндогенные биологические ритмы, которые свойственны всем живым системам. Именно они определяют периодичность деления клеток, отмеряют время рождения и смерти живых организмов [6]. Модели колебательных систем используются в ферментативном катализе, теории иммунитета [7], в теории трансмембранного ионного переноса, микробиологии и биотехнологии [8−10].
Важнейшим примером физиологического осциллятора, период которого составляет порядка одной секунды, является дыхание. Есть и многие другие колебательные процессы с очень коротким циклом, среди которых специфическая нервная активность головного мозга [11].
Несколько отличный тип колебаний наблюдается при гликолизе. Гликолиз — это процесс, в ходе которого происходит расщепление глюкозы для снабжения клетки необходимой для ее жизни энергиейколебания с периодом в несколько минут наблюдаются при определенных концентрациях участвующих в процессе химических веществ [12]. Эти колебания отличны от биологических часов, которые, как отмечено выше, связаны с циркадными или суточными ритмами или внешней периодичностью, в связи с чем их более корректное название — автономные осцилляторы [13]. Колебания в гликолизе и других метаболических системах, периодические процессы фотосинтеза, колебания концентрации кальция в клетке, колебания численности животных в популяциях и сообществах относятся к классу автоколебательных систем. Система называется автоколебательной, если колебания в системе имеют постоянные период и амплитуду, устанавливаются независимо от начальных условий и поддерживаются благодаря свойствам самой системы, а не вследствие воздействия периодической силы [14].
Известным явлением является генерация электрического импульса отдельно взятыми нервными клетками и нейронами. Этому свойству мембран нервных клеток посвящена работа Ходжкина и Хаксли, выполненная на аксоне гигантского кальмара [15]. Анализ явления пространственного распространения потенциалов действия нервных импульсов вдоль по аксону нервной клетки проводится на основе модели Ходжкина и Хаксли: а / а ^ сИ.
С^—8мат" к{У-УМа) + 8кпУ-Ук) + § 1{У-У1) + 1а, где аирзаданные функции потенциала V, g — постоянные проводимости, УМа, Ук, V, — постоянные равновесные потенциалы, С — ёмкость.
Упрощенно сердечную ткань можно рассматривать как среду, состоящую из автоколебательных и возбудимых элементов-клеток. Для моделирования одной сердечной клетки используется биологически релевантная модель Луо-Руди [16] мембранного потенциала кардиомиоцита. Данная модель является моделью типа Ходжкина-Хаксли и описывает изменение мембранного потенциала клетки в зависимости от ионных токов, протекающих через нее [17].
Среди известных периодических реакций одной из важных и базовых является реакция окисления лимонной кислоты броматом калия, катализируемая ионами церия (IV), называемая реакцией Белоусова-Жаботинского [18]. На рис. 1 приведены экспериментальные кривые концентраций реагентов из работы Филда, Кёреса и Нойеса [19]. Сравнительно детальное описание реакции и метод её моделирования будут даны в главе 2 диссертации.
10'.
10 Г.
СеЯ]/[СеШ].
VWW о 300 600 9 00 1200.
Рис. 1. Потенциометрические экспериментальные кривые [Вг ] и 1%([Се1¥-] / [Се III]) для типичных временных осцилляций в реакции Белоусова-Жаботинского.
Практически все модели осцилляторов, описывающих рассмотренные выше явления и процессы, приводят к системам обыкновенных дифференциальных уравнений (ОДУ) для вектора концентраций У (/), а именно, dt F (Y),.
1) где ^(7) = (/1(7),/2(Г),., Л (Г)), У = (>>! (/), у2 (0> Уы (0) • Вектор Р описывает нелинейную реакционную кинетику или механизм, лежащий в основе химических и биологических колебаний [20]. При задании начальных значений система (1) приводит к задаче Копій ОДУ, в общем случае имеющей вид йУ.
— F (t, Y), a t.
Y (t0) = Y0,.
2) где F (t, Y) = (/, (t, Y), f2 (t, Y),., fN (t, Y)), Y = (y, (0, (t), yN (0),.
Ooi> Уо2> ¦¦¦> У on) — При компьютерном приближении решения задачи (2) широко применяются численные методы, которые должны обладать высокой точностью с одновременной минимизацией временной сложности.
Современное состояние вычислительной техники характеризуется непрерывным возрастанием мощности компьютеров по быстродействию и памяти, тем не менее, одним из наиболее важных требований к численным методам, в основном, является минимизация числа операций [21 — 23]. Это связано с необходимостью многократного повторения вычислений, присущей многим алгоритмам решения таких задач, как, например, оптимизация, идентификация параметров, наведение на цели и других [24, 25].
Разнообразие численных методов, предназначенных для решения задачи Копій для ОДУ, показывает, что часто практические требования к методу, среди которых основными являются противоречащие друг другу требования точности и быстродействия, отвергают выбор универсального метода [26, 27]. Это обуславливается, во-первых, тем обстоятельством, что универсализация алгоритма неизбежно ведёт к его усложнению и, соответственно, к увеличению требуемых ресурсов и, во-вторых, из-за высокой вероятности при решении больших по размерности и сложных по построению алгоритмов численного интегрирования задач возникновения проблем, характер которых зачастую заранее предсказать невозможно [28].
Проблема решения задачи (2) сопровождается таким широким спектром изученных не в полной мере особенностей, что теория численных методов решения систем ОДУ выделилась в отдельный раздел вычислительной математики, постоянно пополняющийся новыми результатами. При этом математические утверждения и выводы о предлагаемых методах могут не согласовываться с получаемыми на вычислительных машинах результатами их применения [28, 29]. В связи с этим не меньшее значение при приближении реальных моделей, в частности, описывающих биологические и биохимические осцилляторы, имеет программная реализация методов. Численное интегрирование систем ОДУ во многих современных математических моделях связано с решением таких проблем, как жёсткость и неустойчивость к возмущениям входных параметров [30]. Наличие множества численных методов, с различной эффективностью преодолевающих указанные проблемы (например, методы С. С. Артемьева, Г. В. Демидова, И. Д. Жонголовича, Ю. В. Ракитского, Р. Буллирша, Дж. Стойера, X. Розенброка, Э. Хайрера, Г. Ваннера, Ч. Гира и других авторов), не снижает актуальность создания эффективных, обладающих достаточной простотой программной реализации алгоритмов. В этой связи создание метода, ориентированного, в частности, на приближение задач, моделирующих биологические и биохимические осцилляторы с высокой точностью и одновременной минимизацией временной сложности, а также применимого при решении других классов задач, имеет актуальное значение.
Численные методы для ОДУ. Численные методы для ОДУ естественным образом разделяются на два класса. В один из них входят методы, использующие одно стартовое значение на каждом шаге («одношаговые методы»), а другой образуют методы, опирающиеся на несколько значений решения («многошаговые методы») [29, 31]. Специальные методы разрабатываются для приближения решения жестких систем [28].
Методы Рунге-Кутта. Метод Эйлера для решения задачи Коши для дифференциального уравнения вида [28].
Г = /((, У), У (*о) = Уо (3) а / был описан Эйлером в 1768 году. Метод прост для понимания и программирования: у1+х = у1, у,), для его сходимости не требуются дополнительные, кроме условий существования и единственности решения задачи (3), ограничения на функцию правой части. Рунге (1895) и Хойн (1900) построили новые методы, включив в формулу метода один или два добавочных шага по Эйлеру. Но именно Кутта (1901) сформулировал общую схему того, что теперь называется методом Рунге-Кутта (РК).
Определение. Пусть ^ - целое положительное число («число стадий») и а2> аз1> «32> а, 2> а*,*->—> >-А> с2,.с5 — вещественные коэффициенты. Тогда метод о +сгк, У0+И (а31 к,+а32к2)),.
Ух =Уо + Ь (ь 1 к +.+Ь, кя) называется 5-стадийным явным методом Рунге-Кутта (ЯМРК) для задачи (3). Символическое представление метода (табл. 1) вошло в обычай после статьи Бутчера в 1964 году [31].
Таблица 1.
Табличное представление ЯМРК.
0 агх.
Сз «31 «32 с, О, г ам-1.
А, Ьг Ь,.
В зависимости от выбора коэффициентов и количества стадий различают ЯМРК различных порядков. Хорошие результаты для малого порядка погрешности дают методы Рунге 2, 3 порядка, метод Хойна 3 порядка. Из множества различных формул РК 4-го порядка чаще используется классическая схема РК [28,31]:
Таблица 2.
Классический" метод Рунге-Кутта.
½ 1/2.
½ 0 ½.
1 0 0 1.
1/6 2/6 2/6 1/6.
В некоторых случаях полезно использовать формулы Гилла, Ральстона и Халла, правило 3/8 Кутты и другие методы 4-го порядка [28].
Все описываемые численные методы применимы к приближению решения задачи (2). Сравнение ЯМРК невысоких порядков в [28] производится на примере приближения решения системы.
— = 2 = 10? ехр (5(>, з -1))>, 4, си сИ = О-,) а1 а/ с начальными условиями ^,(0) = 1, / = 1,4. Аналитическое решение системы ух = ехр (зт (72)), у2 = ехр (5зт (72)), у3 = зт (72) + 1, у4= соз (г2) использовано для вычисления абсолютных погрешностей приближения решения системы различными методами. В зависимости от выбора размера шага погрешность приближения уменьшалась до значения 10″ 5:
Рис. 2. Зависимость максимальной глобальной погрешности от числа обращений к подпрограмме вычисления значения функции (/ с).
На рис. 2 сравнивались следующие методы [28]: классический метод Рунге-Кутта (табл. 2), ' ' ' правило 3/8 Кутта, оптимальная формула, и = 0.3587, V = 0.6346, • Ральстон, Халл, и = 0.4, у = 0.45, формула Гилла.
Такие тесты обычно показывают для классической схемы РК слегка худшие результаты, чем для других схем РК, но различия оказываются весьма незначительными. Методы встроены в большинство библиотек программ. В программном пакете МаМСАО ЯМРК четвертого порядка с фиксированным шагом реализует функция г! фхес1 [32, 33], в программном пакете МАТЬАВфункция ос1е23 реализует метод Богацки-Шампайна 3-го порядка [34].
Методы Рунге-Кутта имеют порядок точности р только при достаточной степени гладкости правой части (2), следовательно, повышение порядка метода РК требует наличия соответственного порядку числа непрерывных производных от решения [35]. В большинстве задач, моделирующих колебательные реакции, такую степень гладкости правой части (2) гарантировать невозможно [3, 20].
Для практической оценки погрешности Ричардсон предложил следующую схему. Погрешность значения уг, полученного в результате выполнения двух шагов длины И по ЯМРК порядка р, оценивается по формуле где м> - значение, полученное в результате выполнения одного шага длины 2И. Аппроксимация величины у0+ 2 И) с порядком р +1 может быть выполнена с использованием выражения у2 =у2 +. Идея Ричардсона позволила с одной стороны обеспечить длину шага И, достаточно малую для достижения требуемой точности вычисляемых результатов, а с другой стороныгарантировать достаточно большую длину шага во избежание бесполезной вычислительной работы [36]. Данная схема в отечественной литературе по численному анализу рассматривается как частный случай правила Рунге [37].
Первые вложенные методы РК предложили Мерсон (1957), Ческино (1962) и Зонневельд (1963) [28, 38, 39]. Вместо пользования экстраполяцией Ричардсона они построили такие формулы РК, которые сами содержали бы кроме численного приближенного значения у1 некоторое выражение у1 более высокого порядка. Более экономичные и точные методы такого типа вывели Сарафян (1966), Ингланд (1969) и Фельберг (1968, 1969) [28, 40]. Вложенные методы представляются в следующем виде (табл. 3):
Таблица 3.
Табличное представление вложенных формул РК.
31 «32.
0.1 О, 2.
1 Ь2 Ъ,.
6. ъг Ь,.
При этом величина у{ = у0 + кх±+Ь5 к5) имеет порядок р, а уо+}г (Ь1к1+.+Ь1к1) — порядок который называют порядком «оценщика погрешности» [28]. В табл. 4 приведен один из наиболее распространенных методов 4-го порядка с 6 стадиями (названия вложенных методов приводятся в виде «фамилия р (.
Таблица 4.
Фелъберг 4 (5) О.
¼ 1/4.
3/8 3/32 9/32.
12/13 1932/2197 7296/2197.
1 439/216 -8 3680/513 -845/4104.
½ -8/27 2 -3544/2565 1859/4104 -11/40.
У 25/216 0 1408/2565 2197/4104 -1/5 0.
У1 16/135 0 6656/12 825 28 561/56430 -9/50 2/55.
Метод Дормана-Принса 5 (4) использует несколько иной подход, заключающийся в минимизации членов погрешности для результата старшего порядка [41].
Вложенные формулы с автоматическим управлением величины длины шага получили широкое распространение. В программном пакете МшИСАИ функция rkada. pt, реализующая ЯМРК с адаптацией шага, использует такой механизм управления величиной шага. Данная функция используется для решения медленно меняющихся систем ОДУ, она анализирует скорость изменения решения и соответственно адаптирует размер шага [42]. В МАТЬАВ функция ос1е45 реализует метод Дормана-Принса 5-го порядка [34] Результаты действия механизма управления шагом наглядно демонстрируются при приближении решения системы для модели Лефевера и Николиса [43], называемого «брюсселятор» — базовой модели, являющейся классическим примером автоколебательного поведения концентраций в системе химических реакций.
Предполагается, что шесть веществ участвуют в следующих реакциях:
А—^Х, В + Х—^У + И, 2Х + У-^>ЗХ, из которых вторая является бимолекулярной, а третья — автокаталитической трехмолекулярной реакцией. Согласно закону действующих масс при обозначении концентраций веществ А, В,. через А (х), В (х),. как функций времени х, данные реакции описываются следующими дифференциальными уравнениями:
А' = -кхА, В' = -к2ВХ, О' = к2ВХ, Е' = к4Х,.
X' = к, А-к2ВХ+к3 Х2У-кАХ, Г = к2ВХ-к3Х2У.
Полученная система упрощается: исключаются уравнения О' = к2ВХ и Е' = к4Х, не влияющие на остальные уравнения, предполагается, что концентрации, А и В поддерживаются постоянными и все скорости реакций равны единице. Ввиду выполненных упрощений в обозначениях У (*) = X О), у 2 О) = 70) имеет место система.
ЛУх.
СІХ.
У г.
СІХ, А + уу2-(В + 1) ух, = Ву,-у] у2, которая имеет одну особую точку ^ = -^? = 0 при у.=А, у7=В/А. В йх (1X окрестности этой точки линеаризованное решение неустойчиво только при В>А2+1. Так как все решения ограничены, то при В>А2+1 существует предельный цикл [28] (рис. 3).
1. 2. 3. 4.
Рис. 3. Решения уравнений «брюсселятора», А = 1, В = 3.
Преимущества методов приближения систем ОДУ с адаптацией величины шага особенно проявляются при приближении задачи Коши, а х.
Уг 1 2 «7^ = 3 У1~У1 У 2 а х с начальными условиями у1(0) = 1.01, у2(0) = 3, очень близкими к особой точке [28].
Плотная выдача. Механизм управления длиной шага используемый при численном приближении задачи (2) выбирает точки, в которых вычисляется решение, в соответствии с требованиями к допустимой локальной погрешности. Часто, однако, требуется вычислить и выдать решение в заданных точках, которые к тому же могут быть расположены достаточно плотно. Одна из возможностей получения решения в таких заданных точках состоит в уменьшении длины шага всякий раз, когда это нужно. Такое повторное уменьшение длины шага от почти оптимальной величины нарушает управление длиной шага и может привести к очень сильному росту времени счета и погрешностей округления [44]. С целью обеспечения «плотной» [24] выдачи выводятся непрерывные методы РК, которые способны выдавать значения решения во всех промежуточных точках х' = х0 + в Ь., где 0 < 0 < 1.
Непрерывные расширения методов РК построили Хорн (1983) и Энрайт, Джексон, Нёрсетт и Томсен (1985) [28, 45]. Неплохим инструментом для плотной выдачи результатов и графического представления решения является непрерывное расширение 4-го порядка для метода Дормана и Принса. Такое решение, переходящее в решение ух пятого порядка точности при 0 = 1, определяется следующими формулами [28]:
Ъх (0) =9(1+ 0(-1337/48О + в (1О39/36О + е (-1163/1152)))),.
62(Э) = 0,.
Ь3 (0) = 1ОО02 (1О54/9275 + 0(-4682/27 825 + 0(379/5565)))/3, 64 (0) = -502 (27/40 + 0 (-9/5+ 0(88/96)))/2, ?5 (9) = 18 225 02 (-3/250+ 0(22/375 + 0(-37/600)))/848, 66(0) = -22 02 (-3/10+ 0(29/30+ 0(-17/24)))/7,.
И’о+е'О^о + лЕме)*,.
При высоких ограничениях на точность используют ЯМРК высоких порядков [46]. «Наивысший порядок, фактически достигнутый для явно построенных методов, равен десяти (книга рекордов Гиннеса, с. 333)» [28]. Первые свободные от недостатков вложенные формулы Рунге-Кутта построил Вернер (1978). Как отмечено в [28] прекрасные численные результаты дает метод порядков 8 (7), который предложили Принс и Дорман (1981) [47].
Среди многошаговых методов наиболее употребительны методы интегрирования на сетке с постоянным шагом к при помощи соотношений вида [37].
4).
0 1 = 0 которые принято называть конечно-разностными методами. В вычислительной практике применяют формулы вида (4) со значениями а0 * 0, Ь0 = 0 экстраполяционные и формулы с о0фО, Ь0фО — интерполяционные. В настоящее время из конечно-разностных методов на практике употребляются в осовном методы Адамса [37].
При очень малых значениях допустимой погрешности предпочтительнее применять экстраполяционные программы, благодаря автоматическому увеличению порядка точности. В программном пакете МшкСАИ на методе Булирша-Штера с использованием рациональной экстраполяции основана функция ЪиШоег, которая используется при приближении нежестких систем ОДУ [32, 33]. Обязательное условие использования этой функциидостаточная степень гладкости правой части системы (2).
Из приведенного выше описания следует, что к выбору численного метода определенные требования предъявляет и тип решаемой задачи. Задачи Коши для ОДУ можно условно разделить на мягкие, жесткие, плохо обусловленные и быстро осциллирующие [48]. Специально жестким системам посвящены [49, 50]. К жестким относят задачи химической кинетики, нестационарные процессы в сложных радиоцепях, системы, возникающие при решении уравнений теплопроводности и диффузии методом прямых, и многие другие. Строгого определения понятия жесткости нет, обычно [48] под ним подразумевают наличие в одной системе как быстрозатухающих, так и медленно меняющихся компонент решения. Начиная с 50-х годов для жестких систем стали создавать специальные неявные методы, подробный обзор этих методов дан в [49]. Среди безитерационных, наиболее популярным является метод, предложенный Розенброком [51], который используется при приближении нелинейных автономных систем. Формулы перехода от точки t к 7 = t + т для автономной системы.
5) а t имеют вид S у=у + т^Ь, к&bdquoi=i.
Е-т y"^-)k,=F (y + x tjkj) + х ^ ^ у tJk}, i = l, s,.
U j-i OK j= д F где E — единичная матрица,—матрица Якоби для системы (5), 5 — число ди стадий схемы, аи, уи, Ь, — коэффициенты. Методы Розенброка принадлежат к большому классу методов, которые стараются избежать нелинейных систем и заменяют их последовательностью линейных систем. В [49] такие методы называют линейно неявными методами Рунге-Кутта. В другой литературе такие методы часто называются «полуявными», «обобщенными», «модифицированными», «адаптивными» или «аддитивными» методами Рунге-Кутты [48, 49]. Широко распространены программы ROS4, RODAS, реализующие методы Розенброка четвертого порядка.
Среди итерационных схем лучшими являются неявные 5-стадийные методы Рунге-Кутта. Для автономной системы (5) они имеют следующий вид.
Д S у =у + т^ь, к" к t = F (у + та tJ к j), i = l, s. ,=1 7=1.
Если матрица коэффициентов поддиагональна, то схема явная и не пригодна для жестких задач [48], если суммирование идет до /, то схема называется диагонально-неявной (DIRK). Такие схемы для нахождения каждого к, требуют решения нелинейной системы уравнений, размерность которого равна размерности вектора у. При суммировании до s для нахождения к, приходится решать систему в 5 раз большей размерности, поэтому на практике ограничиваются D/Л/^-методами. В [49] приводятся коды программ RADAU5, SDIRK4, реализующих неявные методы Рунге-Кутта. Также для приближения жестких задач применяются экстраполяционные методы. Широко распространены программы SEULEX и SODEX. В программном пакете MathCAD для приближения жестких задач реализованы функции stiffr, основанная на методе Розенброка, и stiffb, основанная на методе Булирша-Штера. Для их использования предварительно необходимо определить матрицу Якоби для решаемой системы. Такие специальные методы используются, например, при решении задач химической кинетики.
Химические и биологические осцилляторы. Как уже отмечалось выше, колебательные движения являются одной из наиболее быстро развивающихся областей в теоретической биологии [20]. Существует множество полезных математических моделей пространственно однородных колебаний во времени. Считается, что многие клеточные процессы носят колебательный характер и присущее им ритмичное поведение обеспечивает устойчивую основу динамической самоорганизации развития клетки [11]. Последние годы отмечены повышением интереса к исследованию биологических и биохимических систем, в которых появляются незатухающие колебаниявременные и пространственные. Обзорные статьи по колебательным химическим и биохимическим реакциям опубликовали Николис и Портнов, Нойес и Филд [19], Грей и Гольдбер и Каплан [13, 20, 52, 53].
История изучения колебательных реакций началась с работы Лотки (1910), который на примере теоретической реакции впервые показал возможность существования затухающих колебаний. Позже Лотка предложил механизм реакции, носящий теперь имя Лотки-Вольтерра [54]. Следующим историческим событием в области колебательных реакций стала колебательная реакция Белоусова-Жаботинского [18]. Белоусов обнаружил колебания в концентрации катализатора, в его реакции окисления лимонной кислотой броматом в этом качестве выступал церий. Колебания проявлялись в изменении цвета церия, связанными с переходом Се+Ъ в Се+4. Этот эффект ярче проявляется с атомом железа: цвет изменяется с кирпично-красного, когда железо находится в состоянии Fe2е, на ярко-голубой, когда железо находится в состоянии Fe3e. Изучение реакции было продолжено Жаботинским (1964) и сейчас эта реакция известна как реакция Белоусова-Жаботинского. В случае, когда компоненты реакции могут диффундировать, образуются разнообразные сложные структуры. Это явление получило широкий резонанс среди биологов и физиков, интересующихся проблемами пространственно-временной самоорганизации [20] и ее применения в биологическом формообразовании [55]. Сейчас описано немало химических реакций, проявляющих периодическое поведение, и термин «реакция БЖ» применятся, главным образом, к классу реакций, в которых органическое вещество окисляется ионами бромата в присутствии иона металла в кислом растворе [3]. Хотя реакция БЖ представляет скорее химический, чем биологический осциллятор, она рассматривается как прототип биологической колебательной системы.
Модельная система для реакции Белоусова-Жаботинского, предложенная Филд-Кереш-Нойесом [19], количественно воспроизводит реальные химические реакции. Так как теоретические гипотезы могут быть проверены экспериментально, модели реакции БЖ используют как прототип реальной системы. Краткий обзор развернутой реакции БЖ приведен в [56].
Упрощенная модель, предложенная Филдом и Нойесом, получила название «орегонатор» [19]. Схема реакции имеет вид [28]:
А + У-^^Х Х + У— в+х-^>2х+г где А, В — исходные реагенты, Р, (2 — продукты реакции, X, У, 2 -промежуточные соединения: бромистая кислота НВг02, бромид-ион Вг~, Се+А. Концентрации исходных реагентов полагают в модели неизменными.
Уравнение, описывающее изменение концентраций автокатализатора (Х = НВг02), бромид-иона (У = Вг~) и катализатора (г = Се+4) во времени в соответствии с законом действующих масс имеет вид:
7 У г 2 = -к, АУ-к2ХУ + ¡-к, 2, ¿-И 2 5 сИ 3 5.
Численные значения констант скоростей прямых реакций оценены Филдом и Нойесом из экспериментальных данных [19]. Безразмерная форма записи модели «орегонатора» представлена в виде: = з,(г|-г|а + а-да2), йх йх.
Для моделирования [19] используют значения параметров 5 = 77.27, q = 8.375 • 10″ 6, м> = 0.161, / = 1. Модель с данными параметрами: начальные условия ^(0) = ^ у2{0) = 2, >>з (0) = 3, применялась при тестировании в [49] численных методов для решения жестких задач. Система (6) является примером жесткой системы, компоненты решения которой быстро изменяются по величине на много порядков. В [28] отмечено, что «данный пример служит серьезным испытанием для программ численного интегрирования.».
На основании изложенного можно сделать вывод о том, что проблема решения жестких задач с высокой точностью остается актуальной, в частности, в аспекте математического моделирования автоколебательных реакций. Метод, подобранный для аппроксимации с высокой точностью одного типа задач, не встречается при приближении решения задач другого типа. Иными словами, существующие методы приближения решения задачи Коши для ОДУ не инвариантны относительно вида жесткой системы. Необходимо отметить, что существуют проблемы непрерывности полученного приближения, актуальна задача оптимизации распределения узлов интегрирования, в зависимости от гладкости решения на различных участках отрезка интегрирования. Диссертация посвящена уточнению существующих моделей химических и биологических осцилляторов на основе разработки для таких задач специального численного метода и создания комплекса программ, реализующего разработанный метод. При этом требуется, чтобы искомый метод был основан на инвариантном вычислительном алгоритме как для жестких, так и нежестких задач, с целью его компьютерной реализации при = 77.27 0>2+л (1−8.375×10^, ах.
Ф^^оъ-ла+до).
6) = 0.161(^-^3), помощи единых программных процедур. Более точно, формулируется следующая цель.
Целью диссертационной работы является разработка и исследование инвариантного компьютерного метода решения задачи Коши для системы ОДУ с применением к моделированию химических и биологических осцилляторов. Построение метода основано на кусочно-полиномиальном приближении решения в условиях невысокого порядка гладкости правой части при выполнении требования минимизации погрешности и временной сложности, а также требования непрерывности и непрерывной дифференцируемости приближения в случае моделирования процессов с помощью жестких и нежестких систем.
Для достижения поставленной цели в диссертационной работе решаются следующие задачи:
1. Разработать единый компьютерный метод варьируемого кусочно-полиномиального решения задачи Коши для системы ОДУ на основе интерполяционных полиномов Ньютона с числовыми значениями коэффициентов, обеспечивающий высокую точность приближенного решения со свойствами непрерывности и непрерывной дифференцируемости приближения в условиях невысокого порядка гладкости правой части.
2. Доказать сходимость и дать оценку скорости сходимости конструируемого метода в случае его применения для аппроксимации функций и в случае решения с его помощью задачи Коши для системы ОДУ, показать целесообразность его практического применения, выполнить сравнение с методами высоких порядков в аспекте компьютерной реализации математического моделирования процессов с быстро меняющейся динамикой.
3. Выполнить компьютерную реализацию кусочно-полиномиального метода приближенного решения жестких систем ОДУ для уточнения числовых параметров реагентов автоколебательных реакций при моделировании периодической реакции Белоусова-Жаботинского, суточных колебаний концентрации озона в атмосфере и релаксационных автоколебаний в системе гликолиза.
4. Оценить трудоемкость и временную сложность кусочно-полиномиального приближения решения систем ОДУ, указать зависимость трудоемкости метода от точности кусочно-полиномиального приближения решения жестких и нежестких систем ОДУ в процессе моделирования автоколебательных процессов с учетом параллелизма метода.
5. Разработать способ переноса компьютерного метода кусочно-полиномиального приближения решения задачи Коши для системы ОДУ на случай приближенного решения дифференциальных уравнений (ДУ) в частных производных на основе кусочно-полиномиальной аппроксимации функций двух переменных. Выполнить численный эксперимент по точности кусочно-полиномиального приближения и применимости метода для моделирования волновых процессов.
6. Разработать комплекс программ на основе единого метода кусочно-полиномиального решения жестких и нежестких систем ОДУ, адаптируемый к различным классам задач посредством задания числовых параметров подпрограмм, выполнить с его помощью численный и программный эксперимент по сравнению погрешности и временной сложности известных разностных и предложенных кусочно-полиномиальных схем с целью уточнить физические параметры и фазовые портреты математических моделей автоколебательных реакций.
Методы исследования включают вычислительные методы линейной алгебры и математического анализа, численные и аналитические методы решения систем ОДУ, разностные методы решения уравнений в частных производных, методы математического, компьютерного и численного моделирования автоколебательных реакций, элементы теории сложности параллельных алгоритмов, методы теории интерполяции и теории функций вещественной переменной, методы объектного программирования.
Достоверность результатов вытекает из корректного математического обоснования с помощью аналитических оценок погрешности приближений и временной сложности формализованных алгоритмов, подтверждается результатами компьютерного моделирования, программного и численного эксперимента.
Научная новизна результатов диссертационной работы заключается в следующем:
1. Предложен компьютерный метод варьируемого кусочно-полиномиального решения задачи Коши для системы ОДУ, отличающийся от аналогов по построению на основе кусочного приближения решения на подынтервалах интерполяционными полиномами Ньютона с числовыми значениями коэффициентов, а также программной вариацией длин подынтервалов и степеней аппроксимирующих полиномов, что позволяет достигать сравнительно высокой точности при наличии непрерывности и непрерывной дифференцируемости приближенного решения (С. 39 — 47, 53 -64).
2. Показана равномерная сходимость предложенного метода со скоростью геометрической прогрессии к аппроксимируемой функции, а также к решению задачи Коши для системы ОДУ на конечном промежутке из области допустимых значений в условиях двукратной дифференцируемости правой части, что упрощает его применение по сравнению с методами высоких порядков и обеспечивает численное моделирование процессов с быстро меняющейся динамикой (С. 29 — 37, 47 — 53).
3. Реализовано применение кусочно-полиномиального метода для приближенного решения жестких систем ОДУ, на этой основе представлены результаты компьютерного моделирования периодической реакции Белоусова-Жаботинского, суточных колебаний концентрации озона в атмосфере, релаксационных автоколебаний в системе гликолиза. Результаты отличаются повышенной точностью, а также непрерывностью приближенного решения моделирующих систем в допустимых границах трудоемкости, что позволяет уточнить числовые параметры реагентов рассматриваемых автокаталитических реакций (С. 66−91).
4. Выполнены оценки трудоемкости и временной сложности кусочно-полиномиального приближения решения жестких и нежестких систем ОДУ, показана периодическая зависимость трудоемкости от точности кусочно-полиномиального приближения решения при моделировании автоколебательных процессов, а также возможность снижения временной сложности за счет параллелизма метода применительно к компьютерной реализации математических моделей автоколебательных процессов (С. 109 — 117).
5. Разработан способ переноса компьютерного кусочно-полиномиального приближения решения задачи Коши для системы ОДУ на случай приближенного решения ДУ в частных производных. Метод отличается от известных по построению на основе кусочно-полиномиальной аппроксимации функций двух переменных, по точности компьютерного приближения решения линейных гиперболических уравнений, что позволяет его применять для моделирования волновых процессов, описываемых уравнениями данного вида (С. 118−136).
6. Разработан комплекс программ на основе инвариантного метода кусочно-полиномиального решения жестких и нежестких систем ОДУ, отличающийся тем, что адаптация к различным классам задач реализована заданием числовых параметров подпрограмм. При помощи комплекса выполнен численный эксперимент по сравнению погрешности и временной сложности известных разностных и предложенных кусочно-полиномиальных схем с применением узловых интерполяционных значений на основе методов Эйлера, Эйлера-Коши, Рунге-Кутта, Бутчера и Дормана-Принса. Показана меньшая погрешность предложенного метода, что в сочетании с гладкостью приближения позволяет уточнить физические параметры и фазовые портреты математических моделей автоколебательных реакций (С. 39 — 47, 53 — 64, 72 — 77, 154 — 211).
Основные положения, выносимые на защиту:
1. Компьютерный метод кусочно-полиномиального решения задачи Коши для системы ОДУ на основе кусочно-полиномиального приближения решения интерполяционными полиномами Ньютона с числовыми значениями коэффициентов и программной вариацией длин подынтервалов, а также степеней аппроксимирующих полиномов для обеспечения высокой точности приближенного решения со свойствами непрерывности и непрерывной дифференцируемости.
2. Обоснование равномерной сходимости кусочно-полиномиального метода со скоростью геометрической прогрессии к решению задачи Коши для системы ОДУ на конечном промежутке при двукратной дифференцируемости функций правой части с целью обеспечения компьютерного расчета математических моделей периодических реакций.
3. Компьютерная реализация предложенного метода для приближенного решения жестких систем ОДУ, уточнение на этой основе расчета моделей химических автоколебательных реакций, включая числовые параметры реагентов.
4. Оценки трудоемкости и временной сложности кусочно-полиномиального приближения решения жестких и нежестких систем ОДУ для моделей автоколебательных процессовпоказана возможность сокращения времени моделирования за счет параллелизма метода.
5. Метод приближения решения задачи Коши для ДУ в частных производных на основе кусочно-полиномиальной аппроксимации функций двух переменных с применением к расчету моделей волновых процессов.
6. Программный комплекс на основе единого метода кусочно-полиномиального решения жестких и нежестких систем ОДУ, который отличается тем, что адаптация к различным классам задач реализована в виде числовых параметров подпрограмм инвариантного видапри помощи комплекса выполнен расчет математических моделей периодических реакций, получены сравнительные оценки погрешности и времени расчета моделей на основе разностных и кусочно-полиномиальных схем с применением методов Эйлера, Эйлера-Коши, Рунге-Кутта, Бутчера и Дормана-Принсаданы уточнения физических параметров и фазовых портретов математических моделей автоколебательных реакций.
Практическая ценность диссертационного исследования заключается в прикладном характере предложенных методов кусочно-полиномиального решения ОДУ и уравнений в частных производных, которые применяются для компьютерной реализации математического моделирования периодических реакций, включая реакции Белоусова-Жаботинского, суточные колебания концентрации озона в атмосфере, релаксационные автоколебания в системе гликолиза, а также для моделирования волновых процессов. Результаты моделирования необходимы для отладки технологических процессов на основе периодических реакций, для оценки изменений концентрации свободного кислорода, озона и молекулярного кислорода в атмосфере в зависимости от участка земной поверхности. Предложенный кусочно-полиномиальный метод доведен до практической реализации в виде программного комплекса для решения актуальных задач математического моделирования, связанных с исследованием автоколебаний. Кроме того, разработанный программный комплекс на той же основе применяется для снижения временной сложности и повышения точности решения систем ОДУ, моделирующих движение объектов в реальном времени.
Внедрение и использование результатов работы. Полученные в работе результаты использованы:
1.В ОАО НКБ ВС для решения систем ОДУ при моделировании движения транспортного средства в реальном времени (с учетом сил трения и переменного вектора тяги) в трехмерном пространстве. Модель интегрирована в состав программного обеспечения стенда функционального контроля (СФК).
2. В работе по выполнению государственного задания Министерства образования и науки РФ ФГБОУ ВПО «ТГПИ имени А.П. Чехова» по проекту № 7.1398.2011 «Распараллеливаемые компьютерные методы вычисления функций, решения и анализа устойчивости дифференциальных уравнений, цифровой обработки сигналов и распознавания изображений с применением алгоритмов сортировки».
3. В учебном процессе кафедры информатики ФГБОУ ВПО «ТГПИ имени А.П. Чехова» в курсах «Численные методы», «Программирование», «Методы численного анализа и вычислительной алгебры», «Математическое моделирование» и «Компьютерное моделирование».
Апробация работы. Основные результаты работы были представлены на:
— Пятьдесят второй научной студенческой конференции (Таганрог, ТГПИ, 2009 г.);
— III Всероссийской студенческой научно-технической конференции «Прикладная информатика и математическое моделирование» (Москва, МГУП, 2009 г.);
— International Conference Parallel Computer Algebra '2010 (Tambov, Tambov State University named after G.R. Derzhavin, 2010);
— XI международной научно-практической конференции «Фундаментальные и прикладные исследования, разработка и применение высоких технологий в промышленности» (Санкт-Петербург, СПБПУ, 2011);
— XII Всероссийском симпозиуме по прикладной и промышленной математике (весенняя сессия) (Казань, 2011);
— Всероссийской НТК с международным участием: «Компьютерные и информационные технологии в науке, инженерии и управлении» «КомТех-2011» (Таганрог, ТТИ ЮФУ, 2011 г.);
— IX региональной научно-практической конференции «Аспекты развития науки, образования и модернизации промышленности» (Таганрог, ДГТУ, 2011 г.).
Публикации. По материалам работы опубликовано 14 печатных работ общим объемом около 17 печатных листов, в том числе 4 статьи в журналах из перечня рекомендуемых ВАК РФ.
Структура и объём работы. Диссертационная работа состоит из введения, четырех глав основного раздела, заключения, списка литературы и приложений к четырем главам. Основное содержание работы изложено на 152 страницах, включая список литературы из 111 наименований, приложение изложено на 70 страницах, включает коды программ, реализующих математические модели и предложенные численные методы.
4.5. Выводы.
1. В главе на примерах линейных уравнений гиперболического типа показана возможность переноса компьютерного метода кусочно-полиномиального приближения решения задачи Коши для системы ОДУ на случай приближенного решения ДУ в частных производных на основе кусочно-полиномиальной аппроксимации функций двух переменных. Метод отличается от известных по аналитическому и алгоритмическому построению, по точности приближения решения для уравнений рассматриваемого вида, что позволяет рассматривать его применение для численного моделирования волновых процессов, описываемых уравнениями данного вида.
2. Построена численная схема и выполнена программная реализация кусочно-полиномиального приближения решения первой краевой задачи линейных ДУ гиперболического типа, при этом основой для получения узловых значений интерполяционного полинома Ньютона от двух переменных являются сеточные значения разностного решения ДУ в частных производных. Предложенная схема отличается от известных, в частности, тем, что в силу интерполяции и по построению дает непрерывное приближение решения в заданной прямоугольной области.
3. Представлены результаты численного эксперимента по моделированию вынужденных колебаний конечной струны на основе кусочно-полиномиального приближения, которые показывают возможность повышения точности разностных схем в малых по размеру областях.
4. Предложена программная реализация кусочно-полиномиальной схемы приближенного решения линейных ДУ в частных производных гиперболического типа в виде набора стандартных подпрограмм, которые представляют собой основу программного комплекса для численного моделирования волновых процессов, моделируемых уравнениями данного класса.
ЗАКЛЮЧЕНИЕ
.
Основной результат диссертационной работы заключается в разработке и исследовании компьютерного метода варьируемого кусочно-полиномиального решения задачи Коши для жестких и нежестких систем ОДУ с применением к математическому моделированию автоколебательных реакций для повышения точности решения и уточнения концентраций реагентов, а также временных параметров реакций.
В частности, следующие результаты отличаются новизной:
1. Предложен компьютерный метод варьируемого кусочно-полиномиального решения задачи Коши для системы ОДУ, отличающийся от аналогов по построению на основе кусочного приближения решения на подынтервалах интерполяционными полиномами Ньютона с числовыми значениями коэффициентов, а также программной вариацией длин подынтервалов и степеней аппроксимирующих полиномов, что позволяет достигать сравнительно высокой точности при наличии непрерывности и непрерывной дифференцируемости приближенного решения (С. 39 — 47, 53 -64).
2. Реализовано применение кусочно-полиномиального метода для приближенного решения жестких систем ОДУ, на этой основе представлены результаты компьютерного моделирования периодической реакции Бел о у со ва-Жабот ин с ко го, суточных колебаний концентрации озона в атмосфере, релаксационных автоколебаний в системе гликолиза. Результаты отличаются повышенной точностью, а также непрерывностью приближенного решения моделирующих систем в допустимых границах трудоемкости, что позволяет уточнить числовые параметры реагентов рассматриваемых автокаталитических реакций (С. 66−91).
3. Показана периодическая зависимость трудоемкости от точности кусочно-полиномиального приближения решения при моделировании автоколебательных процессов, а также возможность снижения временной сложности за счет параллелизма метода применительно к компьютерной реализации математических моделей автоколебательных процессов {С. 109−117).
4. Разработан компьютерный метод кусочно-полиномиального приближения решения задачи Коши для ДУ в частных производных, отличающийся по построению на основе кусочно-полиномиальной аппроксимации функций двух переменных, а также по точности компьютерного приближения решения линейных гиперболических уравнений, что целесообразно для моделирования волновых процессов (С. 118 — 136).
5. Разработан комплекс программ на основе единого метода кусочно-полиномиального решения жестких и нежестких систем ОДУ, который отличается тем, что адаптация к различным классам задач реализована в виде числовых параметров подпрограмм инвариантного вида. При помощи комплекса выполнен расчет математических моделей периодических реакций, даны сравнительные оценки погрешности и временной сложности расчета моделей на основе разностных и кусочно-полиномиальных схем с применением методов Эйлера, Эйлера-Коши, Рунге-Кутта, Бутчера и Дормана-Принса. Показано, что меньшая погрешность предложенного метода в сочетании с гладкостью приближения позволяет уточнить физические параметры и фазовые портреты математических моделей автоколебательных реакций (С. 39−47, 53−64, 72−77, 154−211).
Работа включает следующие научные результаты.
1. Компьютерный метод варьируемого кусочно-полиномиального решения задачи Коши для системы ОДУ на основе кусочно-полиномиального приближения решения на подынтервалах интерполяционными полиномами Ньютона с числовыми значениями коэффициентов, реализующий программную вариацию длин подынтервалов и степеней аппроксимирующих полиномов для обеспечения высокой точности приближенного решения со свойствами непрерывности и непрерывной дифференцируемости.
2. Обоснование равномерной сходимости предложенного метода со скоростью геометрической прогрессии к аппроксимируемой функции, а также к решению задачи Коши для системы ОДУ на конечном промежутке в условиях двукратной дифференцируемости правой части для обеспечения компьютерной реализации математического моделирования процессов с быстро меняющейся динамикой.
3. Компьютерная реализация кусочно-полиномиального метода для приближенного решения жестких систем ОДУ, уточнение на этой основе результатов компьютерного моделирования различных химических автоколебательных реакций, включая числовые параметры реагентов, за счет повышенной точности и непрерывности приближенного решения моделирующих систем.
4. Оценки трудоемкости и временной сложности кусочно-полиномиального приближения решения жестких и нежестких систем ОДУ при моделировании автоколебательных процессов, а также возможность снижения временной сложности за счет параллелизма метода применительно к моделированию автокаталитических реакций.
5. Видоизменение компьютерного метода кусочно-полиномиального приближения решения задачи Коши для систем ОДУ на случай приближенного решения ДУ в частных производных на основе кусочно-полиномиальной аппроксимации функций двух переменных с применением для моделирования волновых процессов.
6. Программный комплекс кусочно-полиномиального решения жестких и нежестких систем ОДУ, в котором адаптация к различным классам задач реализована заданием числовых параметров подпрограммчисленный и программный эксперимент на основе комплекса по сравнению погрешности и временной сложности известных и предложенных схем с вычислением узловых интерполяционных значений при помощи методов Эйлера, Эйлера-Коши, Рунге-Кутта, Бутчера и Дормана-Принса с целью уточнения физических параметров и фазовых портретов математических моделей автоколебательных реакций.
Научная новизна результатов диссертационной работы заключается в следующем.
1. Предложен компьютерный метод варьируемого кусочно-полиномиального решения задачи Коши для системы ОДУ, отличающийся от аналогов по построению на основе кусочного приближения решения на подынтервалах интерполяционными полиномами Ньютона с числовыми значениями коэффициентов, а также программной вариацией длин подынтервалов и степеней аппроксимирующих полиномов, что позволяет достигать сравнительно высокой точности при наличии непрерывности и непрерывной дифференцируемости приближенного решения (С. 39 — 47, 53 -64).
2. Показана равномерная сходимость предложенного метода со скоростью геометрической прогрессии к аппроксимируемой функции, а также к решению задачи Коши для системы ОДУ на конечном промежутке из области допустимых значений в условиях двукратной дифференцируемости правой части, что упрощает его применение по сравнению с методами высоких порядков и обеспечивает численное моделирование процессов с быстро меняющейся динамикой (С. 29 — 37, 47 — 53).
3. Реализовано применение кусочно-полиномиального метода для приближенного решения жестких систем ОДУ, на этой основе представлены результаты компьютерного моделирования периодической реакции Белоусова-Жаботинского, суточных колебаний концентрации озона в атмосфере, релаксационных автоколебаний в системе гликолиза. Результаты отличаются повышенной точностью, а также непрерывностью приближенного решения моделирующих систем в допустимых границах трудоемкости, что позволяет уточнить числовые параметры реагентов рассматриваемых автокаталитических реакций (С. 66−91).
4. Выполнены оценки трудоемкости и временной сложности кусочно-полиномиального приближения решения жестких и нежестких систем ОДУ, показана периодическая зависимость трудоемкости от точности кусочно-полиномиального приближения решения при моделировании автоколебательных процессов, а также возможность снижения временной сложности за счет параллелизма метода применительно к компьютерной реализации математических моделей автоколебательных процессов (С. 109 — 117).
5. Разработан способ переноса компьютерного кусочно-полиномиального приближения решения задачи Коши для системы ОДУ на случай приближенного решения ДУ в частных производных. Метод отличается от известных по построению на основе кусочно-полиномиальной аппроксимации функций двух переменных, по точности компьютерного приближения решения линейных гиперболических уравнений, что позволяет его применять для моделирования волновых процессов, описываемых уравнениями данного вида (С. 118−136).
6. Разработан комплекс программ на основе инвариантного метода кусочно-полиномиального решения жестких и нежестких систем ОДУ, отличающийся тем, что адаптация к различным классам задач реализована заданием числовых параметров подпрограмм. При помощи комплекса выполнен численный эксперимент по сравнению погрешности и временной сложности известных разностных и предложенных кусочно-полиномиальных схем с применением узловых интерполяционных значений на основе методов Эйлера, Эйлера-Коши, Рунге-Кутта, Бутчера и Дормана-Принса. Показана меньшая погрешность предложенного метода, что в сочетании с гладкостью приближения позволяет уточнить физические параметры и фазовые портреты математических моделей автоколебательных реакций (С. 39−47, 53−64, 72−77, 154−211).
Практическая ценность диссертационного исследования заключается в прикладном характере предложенных методов кусочно-полиномиального решения ОДУ и уравнений в частных производных, которые применяются для компьютерной реализации математического моделирования периодических реакций, включая реакции Белоусова-Жаботинского, суточные колебания концентрации озона в атмосфере, релаксационные автоколебания в системе гликолиза, а также для моделирования волновых процессов. Результаты моделирования необходимы для отладки технологических процессов на основе периодических реакций, для оценки изменений концентрации свободного кислорода, озона и молекулярного кислорода в атмосфере в зависимости от участка земной поверхности. Предложенный кусочно-полиномиальный метод доведен до практической реализации в виде программного комплекса для решения актуальных задач математического моделирования, связанных с исследованием автоколебаний. Кроме того, разработанный программный комплекс на той же основе применяется для снижения временной сложности и повышения точности решения систем ОДУ, моделирующих движение объектов в реальном времени.
Практическое использование результатов работы:
1.В ОАО НКБ ВС для решения систем ОДУ при моделировании движения транспортного средства в реальном времени (с учетом сил трения и переменного вектора тяги) в трехмерном пространстве. Модель интегрирована в состав программного обеспечения стенда функционального контроля (СФК).
2. В работе по выполнению государственного задания Министерства образования и науки РФ ФГБОУ ВПО «ТГПИ имени А.П. Чехова» по проекту № 7.1398.2011 «Распараллеливаемые компьютерные методы вычисления функций, решения и анализа устойчивости дифференциальных уравнений, цифровой обработки сигналов и распознавания изображений с применением алгоритмов сортировки».
3. В учебном процессе кафедры информатики ФГБОУ ВПО «ТГПИ имени А.П. Чехова» в курсах «Численные методы», «Программирование», «Методы численного анализа и вычислительной алгебры», «Математическое моделирование» и «Компьютерное моделирование».
