Содержание к диссертации
Введение
Глава 1. Контроль точности и устойчивости одношаговых методов 10
1.1. Основные определения 11
1.2. Контроль точности вычислений 19
1.3. Контроль устойчивости 23
1.4. Реализация методов с контролем устойчивости 27
Глава 2. Алгоритмы с контролем точности вычислений 30
2.1. Методы типа Рунге-Кутта 30
2.2. Схемы второго порядка точности 35
2.3. Схемы третьего порядка точности 37
2.4. Алгоритм интегрирования на основе метода Мерсона 42
2.5. Анализ результатов расчетов 45
Глава 3. Алгоритмы с контролем устойчивости численной схемы 53
3.1. Схема второго порядка точности 54
3.2. Схема третьего порядка точности 59
3.3. Схема четвертого порядка точности 62
3.4. Анализ результатов расчетов 65
Глава 4. Алгоритмы интегрирования переменного порядка и шага 70
4.1. Алгоритм интегрирования на основе двухстадийной схемы 71
4.2. Алгоритм интегрирования на основе трехстадийной схемы 76
4.3. Алгоритм с применением стадий Рунге - Кутта - Мерсона 81
4.4. Анализ результатов расчетов 88
Глава 5. Комплекс программ RK ODE и результаты моделирования практических задач 94
5.1. Комплекс программ RK ODE для решения жестких задач 94
5.2. Численное моделирование реакции Белоусова-Жаботинского, дающей сложный предельный цикл 100
5.3. Численное моделирование проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма 109
Заключение 115
Приложение
- Реализация методов с контролем устойчивости
- Алгоритм интегрирования на основе метода Мерсона
- Схема четвертого порядка точности
- Алгоритм интегрирования на основе трехстадийной схемы
Введение к работе
В химической кинетике, радиоэлектронике и других важных приложениях возникает проблема численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений
y' = f(t,y), y(t0) = y0, t0
Класс задач, описываемых жесткими системами, расширяется, так как учитывается все большее число факторов при построении математических моделей физических процессов. Поэтому возникает необходимость решения жестких задач все более высокой размерности. Это, в свою очередь, повышает требования к вычислительным алгоритмам. Современные методы решения жестких задач, как правило, на каждом шаге требуют обращение матрицы Якоби, что при достаточно большой размерности задачи определяет общие вычислительные затраты. Известные явные методы, в которых матрица Якоби не применяется, в основном не приспособлены для решения задач даже умеренной жесткости по двум причинам. Во-первых, области устойчивости явных методов малы, что приводит к обременительным ограничениям на величину шага интегрирования. Во-вторых, на участке установления решения шаг раскачивается из-за противоречивости требований точности вычислений и устойчивости численной схемы. Поэтому построение новых эффективных явных методов с расширенными областями устойчивости и контролем устойчивости численных схем, а также алгоритмов переменного порядка и шага, является актушЫнойхэдцдаей.выбрать методы, соответствующие классу решаемых задач. Здесь в основу алгоритмов интегрирования положены явные схемы типа Рунге-Кутта, которые можно записать в виде
Уп+1=Уп+к<РЛ*п,Уп,Ь)> (2)
где п — текущая точка интегрирования, h - шаг интегрирования, cpf —
заданная вектор-функция, зависящая от правой части / задачи (1). Они обладают определенными преимуществами перед многошаговыми
формулами. В частности, многошаговые методы приводят к осреднению решения ("срезание экстремумов"), что при моделировании некоторых динамических объектов делает их неприемлемыми. Обзор работ, посвященных численному решению задачи (1) многошаговыми методами, содержится в [19, 41, 44, 46, 56].
Существует большое количество методов интегрирования жестких систем. Важным і является круг вопросов, связанных с изменением шага интегрирования и оценкой точности получаемых численных результатов, что и делает метод экономичным и надежным. При проведении практических расчетов основным критерием является точность нахождения решения. Поэтому способы управления шагом основаны, как правило, на контроле точности численной схемы. Многие алгоритмы интегрирования при выборе величины шага используют оценку локальной ошибки или, что то же самое, погрешности аппроксимации, потому что если на каждом шаге контролировать некоторый минимальный уровень локальной ошибки, то глобальная ошибка будет ограничена. В настоящее время можно выделить три практических способа оценки данной ошибки ([41], с. 59-65).
Классический способ оценки локальной ошибки одношаговых методов основан на экстраполяционной формуле Ричардсона [81-82]. Его еще называют правилом Рунге. В каждой сеточной точке интервала интегрирования решение вычисляется с шагом h и 0.5/z, а искомая оценка определяется через разность приближений к решению
v0Sh-vh v05h-vh
5" =^ ^ + 0(hp+2), S5Dh=^ ^ + 0(hp+2)
где p - порядок точности метода, yhn и у\ък — приближения к решению с
шагом h и 0.5/z, 8^ и 8^5ph - соответствующие локальные ошибки. Но этот
способ приводит к значительному увеличению вычислительных затрат, т.к. необходимо дважды вычислять решение в каждой точке.
Более дешевым является многошаговый способ [57]. Он заключается в том, что одношаговой формуле /7-го порядка точности в соответствие
ставится многошаговая схема (р +1) —го порядка
/=0
Затем данная формула преобразуется таким образом, чтобы после подстановки в нее приближений (2) получилась оценка локальной ошибки одношагового метода
8п,Р = (X ДГ'Е [«.-^, - Щ/(Уш)1 (3)
;'=0 ;=0
Недостатком этого способа является многошаговость оценки.
В последнее время наиболее популярной является оценка локальной ошибки с помощью вложенных методов. Приближение к решению в каждой точке вычисляется двумя методами вида (2) р—то и (р + 1)-го порядков
точности, а затем локальная ошибка метода р —то порядка оценивается через
разность полученных приближений 5пр = упр+х ~Упр- Такой способ
используется, когда для вычислений по методу р-то порядка не требуется
дополнительных вычислений правой части и матрицы Якоби задачи (1). Следует отметить оперативность и относительную дешевизну оценки локальной ошибки с помощью вложенных методов. По затратам на шаг она лежит между оценкой ошибки с помощью экстраполяции Ричардсона и многошаговой оценкой. В то же время, по отношению к многошаговой оценке, в ней при вычислении ошибки используется информация только с данного шага, что повышает ее надежность. Данный способ успешно применялся в [41, 61, 63, 76, 85-86] и ниже будет использоваться здесь.
Для повышения надежности расчетов необходимо найти оценку глобальной ошибки. Наиболее известный способ определения данной ошибки основан на предположении о линейном характере накопления глобальной ошибки из локальных ошибок [66]. В результате для контроля точности предлагается использовать неравенство
\\SJ
где 5п — оценка локальной ошибки, ||-|| - некоторая норма в RN, є —
требуемая точность расчетов. Такой способ успешно применялся в [56]. Очевидно, что предположение о линейном характере накопления является достаточно грубым. В [3] используется другое неравенство
є/с.
(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 связана с размером области устойчивости. Так как для многих жестких задач длина интервала интегрирования значительно превышает величину D/1 Ятах |, то условие h \ Хтах \< D является весьма и весьма обремКшитселмнБОДов, основанный на дробно — рациональном представлении приближенного решения исследуется в [73], а в [40] рассмотрены методы, основанные на использовании аппарата цепных и ветвящихся дробей. Еще один подход к построению вычислительных алгоритмов заключается в конструировании численных схем, учитывающих специфику исходной задачи. Здесь можно выделить схемы экспоненциальной подгонки [16], и методы, основанные на обращении главной части дифференциального оператора [8-9]. В [15] предложен новый численный метод интегрирования жестких систем, в основе которого лежит принцип последовательной
фильтрации членов, соответствующих наибольшим собственным значениям матрицы Якоби системы (1). Он позволяет без потери устойчивости увеличить шаг интегрирования даже в случае простейших численных схем. В [2] рассмотрены вопросы реализации методов интегрирования на Фортране.
Диссертация посвящена вопросам построения эффективных алгоритмов численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений на основе явных методов типа Рунге-Кутта. Повышение эффективности достигается как за счет более гибкого управления величиной шага интегрирования по точности и устойчивости, так и за счет построения алгоритмов переменного порядка и шага на основе методов с расширенными областями устойчивости. В качестве критерия выбора численной схемы используются неравенства для контроля точности и устойчивости. При решении жестких задач это позволяет на каждом шаге выбирать эффективный с точки зрения вычислительных затрат метод.
Работа состоит из введения, пяти глав, заключения, списка литературы и приложения. Во введении дан обзор работ по теме диссертации и приведено краткое описание содержания диссертации по главам.
Реализация методов с контролем устойчивости
Далее речь пойдет о применении явных методов для решения задач средней или умеренной жесткости. Отметим, что при получении (1.18) может не использоваться информация в точке tn+l, потому что вовсе не обязательно, чтобы даже одна стадия метода вычислялась в точке (tn+l,yn+l). При решении жестких задач в случае быстрого изменения решения это может быть причиной неудовлетворительной точности. Поэтому при получении оценки (1.18) следует привлекать f(tn+liyn+l). В случае успешного шага интегрирования это не будет приводить к существенному увеличению вычислительных затрат, так как f(tn+l,yn+l) используется при выполнении следующего шага. Ниже будем предполагать, что имеется две оценки є п,р и е"п,р глобальной ошибки, причем при вычислении s"n,p используется f(tn+l,yn+l). В принципе, по вычислительным затратам оценка s „tP дешевле чем є"п,Р, но для некоторых методов они могут совпадать.
Теперь для контроля точности вычислений и при выборе величины шага будем контролировать неравенства где - — некоторая норма в RN, є — требуемая точность расчетов. Первое неравенство (1.27) можно использовать для контроля точности вычислений, а при выборе величины шага дополнительно проверять второе неравенство. Более жесткий контроль при выборе шага, чем при контроле точности позволит избежать некоторых повторных вычислений решения (возвратов) вследствие невыполнения точности. В некоторых случаях с целью повышения надежности для контроля точности имеет смысл использовать оба неравенства (1.21).
Обозначим через Vn оценку величины h \ Лтпах , где Хптах есть максимальное собственное число матрицы Якоби. Тогда для контроля устойчивости имеем неравенство Vn D, где D связано с размером области устойчивости метода. Пусть приближение к решению уп в точке tn вычислено с шагом hn. Это означает, что выполнено первое неравенство (1.27), причем, если определить параметр q по формуле qp \\є „,Р\\=є,то q 1. Здесь p — порядок точности метода. Теперь новый шаг hnew{ из условия точности можно вычислить по формуле где параметр s определяется из уравнения sp \\Е"П,Р\\=- При полу&ЬнШ (1.28) использовались условия e nyP = 0{hp) и s"n,p-0{hp). Отметим, что даже в случае успешного вычисления приближения к решению новый шаг интегрирования может быть уменьшен, так как величина s может быть меньше единицы. В этом случае не имеет смысла контролировать устойчивость, а приближение к решению уп+х в точке tn+x вычисляется с шагом hn+l = hnewX. Величину шага hnew2 из условия устойчивости определим по формуле hnew2 = rhn, где параметр гп, учитывая что Vn = 0(К), вычисляется из равенства rVn=D. Теперь, предполагая, что поведение решения не сильно меняется от шага к шагу, величину hn+l можно вычислить по формуле Постоянная 1,1 внесена в (1.29) для того, чтобы создать некоторый зафк гё ) точности и устойчивости и тем самым избежать некоторых повторных вычислений решения. Итак, если неравенства для контроля точности и устойчивости надежны, а спектр матрицы Якоби и производные решения от шага к шагу изменяются незначительно, то приближение к решению уп+х в точке tn+x с шагом hn+l должно быть вычислено успешно. Однако так как оценка максимального собственного числа является грубой, то использование (1.29) приводит к неоправданному колебанию величины шага, что ухудшит эффективность алгоритма интегрирования.
Алгоритм интегрирования на основе метода Мерсона
Одним из самых эффективных среди явных методов типа Рунге-Кутта четвертого порядка точности является метод Мерсона, который имеет вид Пятое вычисление правой части не дает увеличение порядка до пятого, но позволяет расширить интервал устойчивости до 3.5 и оценить величину локальной ошибки 5п4 с помощью kt,l і 5, то есть Во многих библиотеках программ имеется программная реализация метода Мерсона со следующим алгоритмом выбора величины шага интегрирования. Шаг 1. Вычисляются приращения Аг.,1 і 5. Шаг 2. Вычисляется оценка локальной ошибки по формуле (2.39). Шаг 3. Если для некоторого j,l j N, выполняется неравенство 1841 є yJn I, где j есть номер компоненты, то шаг уменьшается в два раза, и снова вычисляются knl i 5. Шаг 4. Вычисляется решение по формуле (2.38). Шаг 5. Если для всех j,l j N, выполняется \S 4\ s\yJn\/32, то шаг интегрирования удваивается. В противном случае остается без изменения. Шаг б: Выполняется следующий шаг интегрирования. В формуле \SJnA\ s\yJn\l32 используется константа 1/32, которая вводится следующим образом. Пусть после успешного шага хотим попытаться увеличить шаг в два раза, то есть hn+x = 2hn. Учитывая, что SnA = 0(h5), можно записать SnA(hn+l) = SnA(2hn) = 25SnA(hn),25=32. Изменение шага интегрирования в два раза обладает некоторыми недостатками. Если при вычислении решения шаг часто и резко меняется, такой алгоритм может приводить к неоправданно большим затратам времени. Кроме того, при неудачном выборе начального шага может понадобиться достаточно много шагов для достижения максимального шага. Обоснование оценки ошибки для контроля точности вычислений проведем на линейном уравнении у = Лу, у(0) = уо, t 0 с комплексным Л,Яе(Л) 0.
Применяя (2.38) для решения данной задачи, получим Уп+\ =0ЗА УП гДе Функция роста Q5A(z) имеет вид 00 z Учитывая, что e: = -, получим, что относительная локальная ошибка схемы (2.38) имеет вид nA-zb/120 + O{z6). Используя формулу (1.20), получим оценку ошибки єп4, то есть Пусть заданная относительная точность вычислений равна є. Учитывая (2.41), имеем z« (720 )025. Теперь можно записать неравенство для контроля точности через оценку локальной ошибки, то есть При получении (2.42) применялось соотношение 5є514 720імє5 4. Несмотря на то, что обоснование (2.42) проведено на линейном скалярном уравнении, оно с достаточно высокой надежностью использовалось для решения нелинейных задач. Кроме того, неравенство типа (2.42) успешно применялось в алгоритмах интегрирования на основе других численных схем [3]. Заметим, что оценка S n,4 = z5/720 + 0(z6) является величиной 0{h5) только в случае f(t,y) = a + bt + cy (см., например, [41], с.62). В других случаях 5п4 содержит некоторые члены порядка 0(h4). Легко убедиться, что для скалярного уравнения (1.7) выполняется соотношение Sn4=-(h4/45)[2f"f3 +6f ff2 +3ff f2] + 0(h5). Заметим также, так как к5 вычисляется в точке tn+l, то отпадает необходимость в построении дополнительного неравенства для выбора величины шага интегрирования. Сформулируем алгоритм переменного шага. Пусть приближение к решению уп в точке tn вычислено с шагом интегрирования hn. Шаг 1. Вычисляется стадия кх по формуле (2.38). Шаг 2.
Вычисляются стадии к2, к2, к4 и к5 по формулам (2.38). Шаг 3. Вычисляется оценка ошибки єп2 по формуле Шаг4. Вычисляется число q\ по формуле q\en4-є, где є есть требуемая точность интегрирования. При определении qx учитывается, что еп4 = 0{hsn). Шаг 5. Если q\ \, то есть требуемая точность не выполняется, то hn полагается равным qxhn и происходит возврат на шаг 2. Шаг 6. Вычисляется приближенное решение уп+] по формулам (2.38). Шаг 7. Вычисляется значение tn+l =tn+h. Шаг 8. Вычисляется шаг по формуле hn+x = qxhn. Шаг 9. Выполняется следующий шаг интегрирования по методу четвертого порядка точности, то есть происходит переход на шаг 1. Ниже алгоритм переменного шага на основе схемы (2.38) будем называть RK4. Отметим, что алгоритм RK4 совпадает с методом Мерсона.
Схема четвертого порядка точности
Как показывают расчеты, контроль устойчивости численной схемы приводит не только к повышению эффективности алгоритма интегрирования, но и позволяет вычислить решение ряда задач, которые не удается рассчитать алгоритмом без контроля устойчивости. Поэтому возникает вопрос о том, как организовать контроль устойчивости известных схем, хорошо зарекомендовавших себя при решении практических задач, не внося существенных изменений в вычислительный процесс. К наиболее известным и удачным явным формулам типа Рунге-Кутта относится пятистадийная схема четвертого порядка точности построенная Мерсоном. Одна из причин успеха данной схемы заключается, по-видимому, в экономичном способе оценки ошибки, на основе которой осуществляется контроль точности вычислений. Кроме того, область устойчивости метода достаточно велика как по вещественной, так и по мнимой оси комплексной плоскости z = Ah. Это позволяет использовать его при решении жестких задач, собственные числа матрицы Якоби которых имеют достаточно большую мнимую часть. На примере этих метода рассмотрим возможность организации контроля устойчивости численных схем с заданными коэффициентами. Для контроля точности метода Мерсона можно применять неравенство Из (3.15) видно, что для данной численной схемы выполняется уо(лЬййф а2 = а3, где а2 = /?21, а3 = /?31 + /?32 Поэтому неравенство для контроля устойчивости получим по аналогии с (3.7). Применяя к разности (к3-к2) формулу Тейлора с остаточным членом в форме Лагранжа первого порядка, будем иметь где вектор уп вычислен в некоторой окрестности решения y(tn). Учитывая, что к2-кх={к21Ъ)/ п/п + 0{кг\ для контроля устойчивости (3.15) можно использовать неравенство где числу 3.5 примерно равна длина интервала устойчивости численной схемы. Отметим, что по мнимой оси область устойчивости (см. рис. 3.3) также ограничена числом 3.5. Сформулируем алгоритм переменного шага с контролем точности и устойчивости численной схемы более детально.
Пусть приближение к решению уп в точке tn вычислено с шагом интегрирования hn. Шаг 1. Вычисляется стадия кх по формуле (3.15). Шаг 2. Вычисляются стадии к2, к3, к4 и к5 по формулам (3.15). Шаг 3. Вычисляется оценка ошибки еп1 по формуле Вычисляется число qx по формуле q[snA=s, где є есть требуемая точность интегрирования. При определении q\ учитывается, что єп4 = 0{h5n ). Шаг 6. Если q\ \, то есть требуемая точность не выполняется, то hn полагается равным qxhn и происходит возврат на шаг 2. Шаг 7. Вычисляется приближенное решение уп+1 по формулам (3.15). Шаг 8. Вычисляется значение tn+l =tn+h. Шаг 9. Вычисляется оценка максимального собственного числа vn 4 матрицы Якоби по формуле vn 4 = 6 шах \к г - к 2\ I \к[ к[\. Шаг 10. Вычисляется число q2 по формуле q2vn4 =3.5. При определении q2 учитывается, что vn 4 = 0{hn). Отметим, что число q2 ограничивает размер шага по устойчивости. Шаг 11. Если qi q\, то есть шаг по устойчивости меньше шага по точности, то новый шаг интегрирования hn+l выбирается равным hn. В противном случае /ги+1 вычисляется по формуле hn+l = mm(qx, q2)h„ Шаг 12.
Выполняется следующий шаг интегрирования по методу четвертого порядка точности, то есть происходит переход на шаг 1. Ниже алгоритм интегрирования на основе численной схемы (3.15), в котором для контроля точности и устойчивости используются, соответственно, неравенства (3.16) и (3.17), будем называть RK4ST. Расчеты проводились для тестовых примеров из химической кинетики, которые приведены в приложении. В качестве критерия эффективности алгоритмов выбрано количество вычислений правой части дифференциальной задачи на интервале интегрирования. Время вычисления не приводится, так как для некоторых задач оно настолько мало, что не может служить объективной характеристикой эффективности метода. Численный эксперимент проводился с целью исследования влияния неравенства для контроля устойчивости на эффективность алгоритма интегрирования. Расчеты проводились на Pentium 4 (с двойной точностью).
Результаты расчетов приведены в табл. 3.1 - 3.24. примерно одинаковое. Число возвратов для методов с контролем устойчивости значительно меньше. Отметим, что. возвраты непременно связаны с дополнительными вычислительными затратами. Поэтому затраты, которые выражаются в количестве вычислений функции /, для методов с контролем устойчивости ниже, чем для методов без контроля устойчивости. Глава 4. Алгоритмы интегрирования переменного порядка и шага В данной главе на основе явных двухстадийной, трехстадийной и пятистадийной численных схем типа Рунге-Кутта разработаны алгоритмы переменного порядка и шага для решения задач средней жесткости. Выбор эффективного метода осуществляется на каждом шаге по точности и устойчивости.
Алгоритм интегрирования на основе трехстадийной схемы
Впервые колебательную химическую реакцию, проявляющуюся в виде периодических вспышек при окислении паров фосфора, наблюдал Роберт Бойль в конце XVII века. Эти повторяющиеся вспышки затем неоднократно описывали многие исследователи. В XIX веке были обнаружены и другие колебательные реакции. Однако они не привлекли особого внимания, поскольку химическая кинетика как наука еще не существовала, и никто не имел представления о том, как должна идти химическая реакция. Лишь во второй половине ХГХ века возникли термодинамика и химическая кинетика, положившие начало специфическому интересу к колебательным реакциям и методам их анализа. Предсказания возможности колебаний в химических системах делались, начиная с 1910 года, на основе анализа системы дифференциальных уравнений, можно найти в работах А.Лотки. Однако первые математические модели соответствовали неосуществимым и невозможным химическим реакциям. К тому же все попытки экспериментально обнаружить колебательные реакции долгое время не давали положительных результатов.
В 1921 году У. Брей опубликовал статью, в которой достаточно подробно описана первая колебательная жидкофазная реакция разложения пероксида водорода, катализируемая иод атом. Хотя эксперимент осложнялся выделением кислорода, Брей осознал связь между своим открытием и прогнозом Лотки. Однако его работа не вызывала интереса в течение примерно 40 лет. Одна из причин такого безразличия - довольно низкий уровень развития методов исследования механизмов сложных химических реакций. Другой причиной было широко распространенное мнение, что второй закон термодинамики запрещает такие колебания даже вдали от равновесия. Фактически большинство химиков считали, что колебания концентрации в закрытых гомогенных системах невозможны, иначе говоря, чисто химических колебаний не бывает.
Современная история исследований колебательных химических реакций в жидкой фазе началась в 1951 году, когда Б.П. Белоусов открыл колебания концентраций окисленной и восстановительной форм церия в реакции взаимодействия лимонной кислоты с броматом, катализируемой ионами церия. Раствор регулярно менял свою окраску от бесцветной к желтой, затем снова к бесцветной и т.д. Белоусов провел достаточно подробное исследование этой реакции и, в частности, показал, что период колебаний сильно уменьшается с повышением кислотности среды и температуры. Реакция была удобна для лабораторных исследований.
Колебания можно было легко наблюдать визуально, а их период находился в пределах 10—100с, совпадая с естественным масштабом времени человека наблюдателя. В то время как химики, к которым присоединились и биохимики, дружно отвергали химические колебания, последние продолжали привлекать внимание математиков и физиков, интересовавшихся биологией. В 1952 году появилась статья A.M. Тьюринга "Химические основы морфогенеза", в которой он показал, что сочетание химических колебаний с диффузией молекул может приводить к появлению устойчивых пространственных структур, где области высоких и низких концентраций чередуются.
Со временем становилось ясно, что второй закон термодинамики не нарушается в живых системах и не мешает их сложному поведению и эволюции. Но для существования жизни или любой ее физической или химической модели необходимо, чтобы система достаточно долго находилась вдали от термодинамического равновесия. Гомогенные химические системы оказались удобной моделью для того, чтобы перейти от общих рассуждений к конкретному анализу.
В 1955 году Илья Романович Пригожий, бельгийский физик и физико-химик, лауреат Нобелевской премии 1977 года за работы по термодинамике необратимых процессов, показал, что в открытой системе химические колебания возможны около стационарного состояния, достаточно удаленного от термодинамического равновесия; колеблется только величина скорости изменения энтропии, а ее знак всегда остается положительным.
В работах И.Р. Пригожина и его сотрудников было показано, что колебательные химические реакции не только вероятны, но и . возможны Старая парадигма, утверждающая, что Природа запрещает колебательные реакции, сменилась новой, в которой они рассматриваются как интересная и фундаментально важная область науки. С этих работ начался современный этап исследований колебательных химических реакций.
В конце 1961 года работа Б.П. Белоусова была продолжена A.M. Жаботинским, который получил колебания при использовании в качестве восстановителя в реакции Белоусова не только лимонной, но и малоновой, а также яблочной кислот. A.M. Жаботинский провел подробные исследования колебаний в системе с малоновой кислотой, которая оказалась более удобным восстановителем, так как протекание реакции не осложнялось газовыделением. Новость об этой изумительной реакции обошла весь мир, и в нескольких лабораториях (в СССР, США и Западной Европе) стали интенсивно изучать реакцию, которая сейчас широко известна под названием "реакция Белоусова-Жаботинского". Впоследствии A.M. Жаботинским было показано, что автоколебательная реакция может осуществляться и в том случае, когда лимонная кислота заменена малоновой или любой другой дикарбоновой кислотой с активной метиленовой группировкой, а каталитическая редокс-пара Се(ГУ)/Се(Ш) заменена парой Мп(Ш)/Мп(П) или ферроин / ферриин. Например, при взаимодействии Се4+ с малоновой кислотой в присутствии бромат-иона раствор периодически окрашивается то в вишнево-красный цвет, то становится бесцветным. Такие химические часы будут тикать до тех пор, пока в системе есть бромат и малоновая кислота.
