О повышении точности численного решения дифференциальных уравнений химической кинетики Текст научной статьи по специальности «Математика»

УДК 519.673

В. Х. Федотов, Н. И. Кольцов

О ПОВЫШЕНИИ ТОЧНОСТИ ЧИСЛЕННОГО РЕШЕНИЯ

ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ХИМИЧЕСКОЙ КИНЕТИКИ

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

химическая кинетика.

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

Keywords: ordinary differential equations, Cauchy problem, numerical methods, approximate calculations, chemical kinetics.

New simple numerical methods for solving of ordinary differential equations (ODE) were developed. The efficiency of their application in comparison with known numerical methods on examples solving of ODE chemical kinetics was showed.

Введение

Известно [1], что химическая кинетика - наука о скоростях и закономерностях протекания химических реакций во времени. Один из основных этапов ее исследования связан с решением задачи Коши -численным интегрировании обыкновенных дифференциальных уравнений (ОДУ) с начальными условиями, описывающих зависимости концентраций реагентов во времени. От точности используемого метода решения зависит адекватность описания кинетики конкретной реакции определенными ОДУ. Особенно это важно для реакций, характеризующиеся сложным кинетическим поведением [2-4]. В данной статье разработаны методы повышения точности численного решения задачи Коши без увеличения порядка приближения с иллюстрацией их эффективности по сравнению с аналогичными известными численными методами на примерах решения ОДУ химической кинетики конкретных химических реакций.

Результаты и их обсуждение

Рассмотрим задачу Коши для системы ОДУ

вида

у = с1уМ = ШШ У(Ы = Уо, (1) где у(Т) = (у1, у2,..., у к) - вектор зависимых переменных, у0 - вектор начальных значений переменных, t - независимая переменная (время), ^ - начальный момент времени (начало эволюции), f -известная (заданная) функция. Задача Коши в общем виде не имеет точного решения в квадратурах [5]. Численные методы решения основаны на различных разностных схемах, представляющих собой конечное приближение бесконечного ряда, аппроксимирующего функцию уф. Чем выше порядок приближения, тем, как правило, ниже погрешность метода (увеличение порядка приближения не обязательно снижает погрешность). Целью данной работы является построение простых методов численного решения задачи Коши с минимальной погрешностью без увеличения порядка приближения, которые можно отнести к модификациям методов Эйлера.

1. Одномерные системы. Все модификации методов Эйлера основаны на аппроксимации зависимой переменной у(0 конечным числом p членов бесконечного ряда Тейлора в малой окрестности начального приближения [5]

уЦоЩ * у{Ь) + У № + у"(^)Фг/2\ +

у'" (10)Ф3/3\ +.... + ур)%)Ф /р\................(2)

Погрешность формулы (2) (остаток порядка р) определяется бесконечным рядом

Яр(Ф) = у/р+1){Ь)^Р+1 /(р+1)\ + у(р+2)%)^+2 /(р+2)\ +

... ~ О(р) ~ 0фр+1).............................................(3)

Здесь 0фр+^) обозначает бесконечно-малую порядка р+1 [6], т.е. < Сгде С ~ 1, а о(Фр) - бесконечно-малую более высокого порядка, чем р, т.е. < &Фр+1, где е<<1. Основной вклад в эту погрешность вносит первое слагаемое, зависящее от (р+1)-й производной в точке ^ и шага h. Сумма остальных образует погрешность порядка малости р+2 и т.д.

При разработке алгоритмов численных вычислений аналитические соотношения представляют в виде дискретных итерационных формул. Подставляя в (2) вместо у ее известное значение /из (1) с учетом (3) при р = 1, получим обобщенную разностную формулу первого порядка:

уп+1 = уп + Ф у (П + ЪФ) * уп + Ф гуп, и)+

п=0,1,2,.... , ......................................................(4)

где уп - значение на предыдущей итерации, уп+1 -значение на следующей итерации, - локальная ошибка на текущей итерации, 0<Ф<<1 - шаг разностной схемы. В этом случае основной вклад в погрешность вносит вторая производная в точке tn. При ¡31ф) = 0 из (4) следует классический метод Эйлера первого порядка, для которого локальная ошибка максимальна = Яэйлера ~ 0(Ф ).

На практике, погрешность итерационных формул вида (4) зависит от мало изученных законов алгебры приближенных вычислений (нарушение коммутативности и др.), определяемых локальной ошибкой метода (на данной итерации), ошибками дискретизации (округлений, отсечений), устойчивостью разностных схем и др. [7]. Их следствием является то, что вычисление уг+1 для каждой итерации фактически означает решение задачи Коши, отличной от предыдущей, т. к. представляет собой линию

семейства кривых y(t, Const) с началом в точке yn. Этот эффект «сползания» может привести к накоплению больших глобальных ошибок R(h) = ZRP(h) на всем интервале интегрирования. Например, глобальная ошибка метода Эйлера на порядок выше локальной R(h) ~ O(h) [7]. Бесконечный остаток Rp(h) в соотношении (3) является ошибкой метода, а его конечное приближение R1n в (4) является ошибкой дискретизации. Они близки, но не равны и обе «немного» отличаются от неизвестной истинной (точной) локальной ошибки En решения исходной задачи Коши (1), т.к. неизвестно ее точное решение y(tn). Критерий высокой точности численного решения задачи Коши запишем в виде задачи нелинейной минимизации локальной ошибки на интервале

[tm+1, tn]:

Ri,n * En = I yn+i-y(tn) I ^ min. (5)

Соотношения (4), (5) за счет выбора разных способов оценки остатка позволяют конструировать алгоритмы численного интегрирования различной точности. Например, оценки остатка в форме Лагранжа и Коши используют значение второй производной внутри отрезка (tm+i,tn) интегрирования [8]

Rl,n * Rлагранж = У"(6п)^/2 ~ O(h2), (6)

Ri,n * RKomn = y"(0n)(1-9n)h2 ~ O(h2), (7) где 6n = tn+anh, 0<an<1. Эти оценки используют значения высших производных в неопределенных точках и являются разновидностями теорем о среднем [5], которые недостаточно неконструктивны. Если высшие производные можно найти из условий задачи (1) через полные дифференциалы y" = fyf + ft, y" ' = ff + fy + ftt, то an , 8n , ßn остаются неизвестными.

Найдем боле точные и конструктивные формы оценки R1n. Для этого рассмотрим поведение y(t) вблизи точки tn. Локально, любую непрерывную и дифференциируемую функцию можно считать монотонной. Возможные случаи: 1) y(t) выпукла вниз (лежит над касательной); 2) y(t) выпукла вверх (лежит под касательной); 3) y(t) имеет перегиб (лежит по разные стороны от касательной).

Случай 1. Пусть функция y(t) на интервале [tn, tn+1] лежит над касательной (выпукла вниз, y (tn) ^ 0, y" (tn) > 0). Тогда метод Эйлера на каждой итерации занижает значение yn+1 по сравнению с точным значением y(tn) на неизвестную величину E = 1 yn+1- y(tn)l * R n. Оценим ее через первые производные, значения которых y (tn) = f(tn) для задачи Коши известны в любой предыдущей точке. Грубой оценкой является разность между двумя последовательными приближениями R1n < Iyn+1-ynI = h I y (tn) I. Улучшим оценку за счет неизвестного значения первой производной на правом конце tn+1 интервала R1n< h I y (tn)I< h I y (tn+1)I. Эту оценку можно улучшить за счет другого значения первой производной в точке 6n = tn + an(tm+1-tn) = tn + anh = antn+1 + (1 -an)tn, 0 < an < 1, расположенной немного левее an = 1-е правого конца интервала R1n < h I f(6n) I. Понятно, что улучшение можно продолжать, сдвигаясь еще левее. При этом внутри интервала должно существовать неизвестное оптимальное

значение 6n*, обеспечивающее минимум ошибки En (разновидность теорем о среднем) [5]

h min I f(tn+an*h) I < Яоптим * h I/(9n*)I < h maxi f(tn+a.n*h) I, (8) где 6n* = tn + an*h - оптимальное значение, 0 < an* < 1 - коэффициент оптимальной точности. Оценка (8) в отличие от (6) использует только первую производную и не зависит от производных высших порядков. Чтобы сделать оценку (8) конструктивной, зададим алгоритм вычисления / (9n*), например, с помощью линейной / (6n*) * an*y (tn) + yn интерполяции [5]. Для линейной интерполяции при yn = 0 соотношение (4) с учетом (8) запишется в виде оптимизированного алгоритма Эйлера

Уп+1 * Уп + hyn' + Коптим = Уп + hyn' + an*hyn' =

Уп + h(1 + an*)/n', n=0,1,2,...........................(9)

Критерием точности выберем разность между интерполированным значением и нужным числом членов остатка. Для определения нужного числа членов остатка Rp(h), зададим желаемую точность ~ O(hk). Тогда критерий (5) запишется в виде уравнения En*Ron™M-Rfc(h)*h[an*Iyn' I -Zy(p+1)n h/(p+1)!]=0,

p=1,2,3,...,k-1. (10) корень которого дает оптимальные значения an*:

an* = Zy(p+1)nhp /[(p+1)! I yn I], p=1,2,23.....k-1. (11)

Так, чтобы обеспечить точность ~O(h ), полагаем k

= 2, т.е. p = 1 и из (11) получим an* = y"nh / (21 yn'I).

3

Чтобы обеспечить точность ~O(h ), полагаем k = 3 и из (11) получим an* = [y""nh/2 + / "nh2/6]/1 yn 'I и т.д. Отметим, что критерий точности (5), в отличие от самого метода (9), использует производные высших порядков. Однако, практика показывает, что достаточно применить его один раз на первом (самом грубом) шаге. В численных вычислениях критерием точности обычно считают абсолютную величину разности двух последовательных приближений Дп = I yn+1 —yn I, которая легко вычисляется, но, отличается от истинной погрешности En и может неверно оценивать реальные ошибки вычислений. Критерий (10), в этом смысле, конечно, объективнее.

Возможны и другие способы определения an*. Простейшими методами являются перебор и различные эвристики. Например, для эвристики «среднего» acpe^ = 1/2 и оценка (8) запишется

R1,n * ^среднее = 0.5h I/n 'I. (12) Для эвристики «золотого» сечения 1/a = a/(1— a) существует два симметричных значения a30r]0T * 0.382, 1-a30]0T * 0.618 и два варианта записи оценки (8) Кзолото1 = 0.382h I y nI, Кзолото2 = 0.618hI y 'nI. Сравним полученные результаты с методами Эйлера первого (Эйлер 1) и второго (Эйлер 2) порядков.

Пример 1. Рассмотрим необратимую реакцию А^В, нестационарное поведение которой, согласно закону действующих масс, описывается уравнением y = -ky, (13)

где у - концентрация вещества А, k>0 - константа скорости реакции. Задача Коши для уравнения (13) заключается в расчете значения у на отрезке времени [0,h] при начальных значениях t0 = 0, y0 = 1. Ее

точное решение у(0 = е-к выпукло вниз (лежит над касательной) и зависит от к.

Пусть k = 10. В табл. 1 приведены результаты численного интегрирования уравнения (13) известными (Эйлер 1, Лагранж, Эйлер2) и предлагаемыми нами (Среднее, Золото 1, Оптим) методами. Из строки 2 табл. 1 видно, что наименьшую локальную ошибку на первом шаге дает оптимальная модификация ^оптим = 0.00002 ~ 0(Н5) при значении аоптим ~ 0.3679, найденном бисекцией. Оптимальное значение ап*, определенное по формуле (11) для р=1, равно ап* = -^/2 = 1/2, что совпадает с модификацией среднего. Для р = 2 а/ = -Ш2+ к^Ь2/6 = 1/21/6 = 1/3, для р=3 а* = -Ш2+ №/6- a3h3/24 = 1/3+1/24 = 9/24 и 0.375, для р=4 аn* = -Ш2+ №/6- к3^/24 + к4Ф/120 = 1/3+1/24 = 9/24-1/120 = 44/120 и 0.366667, для р=5 аn* = -Ш2+ к2Н2/6-к3Ф/24+к4Ф/120 - к5Н5/720 = 44/120 + 1/720 = 265/720 = 0.368056, для р=6 аn* = -Ш2+ к2Н2/6-к3Ф/24+ к4Ф/120 - к5Ф/720+ к6Ф/5040 = 265/720-1/5040 = 0.367857 и аоптим. Это значение обеспечивает нужную точность. Ошибка метода «Среднее» Есреднее = 0.4 ~ Оф), который также не использует вторую производную, сравнима с ошибками по Лагранжу и методу Эйлера второго порядка, Наименьшую глобальную ошибку также дал оптимальный метод, даже без пересчета оптимального значения на каждом шаге. Если на каждом шаге пересчитывать аоптим, то результаты будут еще точнее, см. рис. 1. Расчеты значений у(0 известными и предлагаемыми нами методами при разных значениях шага h приведены в табл. 1. Как видно, с уменьшением шага до 0.01 точность расчета возросла. Оптимальное значение аоптим и 0.049 обеспечивает на втором шаге Еоптим = 0.00006 ~ O(h ) и гло-

2

бальную точность £оптим = 0.0001 ~ O(h ).

0,4

:

Эйлер

Эйлер2 •

-Оптим

Рис. 1 - Зависимости у({), полученные различными методами интегрирования для уравнения (13) при к = 10, Л = 0.1, у0 = 1

Таблица 1 - Точность расчета различными методами для уравнения (13) при к = 10

t Точн Э2 А Сред А Зол1 А Опт А

h=0.1

0.0 1.00 1.00 0.00 1.00 0.00 1.00 0.00 1.00 0.000000

0.1 0.37 0.50 0.13 0.50 0.13 0.38 0.01 0.37 0.000021

0.2 0.14 0.25 0.11 0.25 0.11 0.15 0.01 0.14 0.000015

0.3 0.05 0.13 0.08 0.13 0.08 0.06 0.01 0.05 0.000008

0.4 0.02 0.06 0.04 0.06 0.04 0.02 0.00 0.02 0.000004

0.5 0.01 0.03 0.02 0.03 0.02 0.01 0.00 0.01 0.000002

0.6 0.00 0.02 0.01 0.02 0.01 0.00 0.00 0.00 0.000001

Итог 0.40 0.40 0.04 0.000051

h=0.01

0.00 1.00 1.00 0.0000 1.00 0.00 1.00 0.00 1.00 0.000000

0.01 0.90 0.91 0.0002 0.95 0.05 0.94 0.03 0.90 0.000063

0.10 0.37 0.37 0.0007 0.60 0.23 0.53 0.16 0.37 0.000255

0.20 0.14 0.14 0.0005 0.36 0.22 0.28 0.14 0.14 0.000187

0.30 0.05 0.05 0.0003 0.21 0.16 0.15 0.10 0.05 0.000103

0.40 0.02 0.02 0.0001 0.13 0.11 0.08 0.06 0.02 0.000051

0.50 0.01 0.01 0.00006 0.08 0.07 0.04 0.03 0.01 0.000023

0.54 0.00 0.00 0.00004 0.06 0.06 0.03 0.03 0.00 0.000017

Итог 0.0004 0.47 0.23 0.000153

Примечание: Э2 - метод Эйлера второго порядка, Сред - метод среднего, Зол1 - метод золотого сечения, Опт - оптимальный метод, Д - ошибка.

Расчет по методу Эйлер1 дает сильно заниженное значение у = 0.0. Расчет по методу Эйлер2 требует вычисления второй производной у" = у + Ъ = ^у, но дает завышенное значение. Таким образом, предлагаемый нами метод «Оптим» позволяет точнее рассчитать значения у(И) по сравнению с другими методами не зависимо от значения шага интегрирования h.

Случай 2. Если функция уф лежит под касательной (выпукла вверх, у ^п) Ф 0, У"($П) < 0), то все происходит с точностью до «наоборот». На каждой итерации метод Эйлера завышает значение уп+1 по сравнению с точным у^П) на ту же величину, и в алгоритме (9) нужно поменять знак перед ап*

Уп+1 и Уп + - ап*)уп, п=0,1,2,..... (14)

Для удобства записи соотношения (9) и (14) можно записать в форме одного общего алгоритма

Уп+1иУп+^1+В\дп(Уп'")ап*)Уп', п=0,1,2,...(15) Заметим, что в этой форме алгоритм использует знак второй производной. Расчеты показали, что и в этом случае наилучшие результаты показывает оптимизированный метод.

Случай 3. Пусть функция у(И) имеет перегиб (лежит по разные стороны от касательной, у ф) Ф 0, У"((п) = 0). Тогда справа и слева от точки перегиба решение лежит над касательной (для четных функций) и под касательной (для нечетных функций). Для прямого времени (эволюция вперед) определяющим является поведение справа (для обратного времени - слева). Для этого вырожденного, но не типичного случая, алгоритм (14) можно скорректировать так

уп+1 иуп+^1+э\дп(у"(уп+Ы2))ап*)уп', п=0,1,2 (16) 2. Многомерные системы. Все приведенные рассуждения обобщаются на многомерный случай, т.е. на системы ОДУ, аналогично классическим многомерным методам Эйлера. При этом скалярное соотношение (15) следует рассматривать в векторной форме

уп+1.1и уп. , + ^1+в\дп(уп, , '")ап*)уп, , ',

п=0,1,2,....; ¡=1,...,к. (17)

Анализ показал, что точность найденных модификаций метода Эйлера (8)-(17) превосходит известные методы, например Рунге-Кутта, имеющего порядок точности ~ 0(Л5). Рассмотрим примеры применения модификации метода Эйлера (17) для многомерных систем.

Пример 2. Рассмотрим реакцию А=В, протекающую по трехстадийной схеме [9] 1) А + Х3 = X!, 2) Х-| = Х2, 3) Х2 + 2Х3 ^ 3Х3 + В, нестационарное поведение которой при постоянстве концентраций основных веществ А и В описывается системой уравнений

х' 1 = ю1(1-х1-х2)-ю_1х1-ю2х1+ю_2х2, X 2 = ю2х1-ю-2х2-ю3х2(1 -х1-х2)2, (18)

где х, - поверхностные концентрации промежуточных веществ Х,; ю1 = к1СА, ю_1 = к-1, ю2 = к2, ю_2 = к-2 и ю3 = к3 - частоты стадий, к+, -константы скоростей стадий, СА - концентрация вещества А. Динамика системы (18) может быть монотонной, колебательной (затухающей и незатухающей), т. е. сочетает участки выпуклости вверх и вниз. Стационарных состояний может быть три или одно. Если стационарное состояние единственно и неустойчиво, то система (18) описывает незатухающие колебания (автоколебания), что реализуется, например, при значениях частот стадий ю1 = 2.89, ю-1 = 0.01, ю2 = 3/89, ю-2 = 0.1, ю3 = 2000, которым соответствует неустойчивое стационарное состояние х1м = 0.89; х2м = 0.1. Результаты численного интегрирования системы (18) с помощью различных модификаций метода Эйлера (Эйлер 3 - третьего порядка, Точно -максимальная точность) и разработанного нами метода (Оптимальный) при начальных значениях х10 = 0.7, х20 = 0.2 и Л = 0.01 приведены в табл. 2.

Таблица 2 - Точность расчета значений х1 (/) различными методами для системы (18)

t Э2 Д Э3 Д Точн Д Опт Д

0.00 0.70 0.70 0.70 0.70

0.01 0.70 0.0000 0.70 0.000 0.70 0.00 0.70 0.0052

0.02 0.71 0.0001 0.71 0.000 0.71 0.00 0.70 0.0049

0.03 0.71 0.0002 0.71 0.000 0.71 0.00 0.71 0.0070

0.10 0.76 0.0003 0.76 0.000 0.76 0.00 0.76 0.0035

0.20 0.82 0.0004 0.82 0.000 0.82 0.00 0.82 0.0026

0.30 0.86 0.0005 0.86 0.000 0.86 0.00 0.86 0.0019

25.00 0.67 0.3075 0.88 0.013 0.88 0.00 0.88 0.0001

25.20 0.67 0.3140 0.88 0.017 0.89 0.00 0.89 0.0002

25.40 0.67 0.3238 0.88 0.024 0.89 0.00 0.89 0.0003

25.49 0.67 0.3298 0.88 0.027 0.89 0.00 0.89 0.0004

Д 0.1907 0.006 0.00 0.0002

= 4ю3(1-х1-х2) - 2со3х2 и х2"" = - 6со3, х2 '"" = 0. Как видно, наилучшие результаты показывает оптимизированный метод, среднеквадратическая ошибка (корень из суммы квадратов, деленный на число точек) минимальна Еоптим = 0.0002 даже при не полностью оптимальном значении а1 = а2 = 0.5. Графические иллюстрации приведены на рис. 2. Как видно из этого рисунка, увеличение порядка приближения до двух (Эйлер2) увеличило погрешность.

1 ,00

□ 5 10 15 20 25

Рис. 2 - Зависимости Xi(f), полученные для системы (18) различными методами

Таким образом, нами разработаны новые простые методы «Среднее», «Золото1», «Золото2» и «Оптим» численного решения задачи Коши для обыкновенных дифференциальных уравнений без использования производных высших порядков. Метод «Оптим» позволяет интегрировать дифференциальные уравнения с максимально высокой точностью и минимальным эффектом «сползания». Он может быть применен для уточнения решений систем ОДУ, в частности, уравнений химической кинетики для более точного расчета траекторий изменения концентраций реагентов реакций, характеризующихся сложным нестационарным поведением.

Литература

1. С. Л. Киперман. Основы химической кинетики в гетерогенном катализе. М., Химия, 1979, 349 с.

2. В.Х. Федотов, Н.И. Кольцов, Вестник Казан. технол. ун-та, 16, 23, 7-9 (2013).

3. В.Х. Федотов, Н.И. Кольцов, Вестник Казан. технол. ун-та, 16, 23, 10-12 (2013).

4. В.Х. Федотов, Н.И. Кольцов, Вестник Казан. технол. ун-та, 16, 24, 17-21 (2013).

5. Г. Корн, Т. Корн. Справочник по математике для научных работников и инженеров. М., Наука, 1974, 832 с.

6. В.П. Дьяконов. Справочник по алгоритмам и программам на языке бейсик для ПЭВМ. М., Наука, 1989, 242 с.

7. Дж. Форсайт, М. Малькольм, К. Моулер. Машинные методы математических вычислений. М., Мир, 1980, 280 с.

8. Г.М. Фихтенгольц. Курс дифференциального и интегрального исчисления, т. 1. М-Л., ОГИЗ Гостехиздат, 1949, 690 с.

9. V.Kh. Fedotov, B.V. Alekseev, N.I. Koltsov, React. Kinet. Catal. Lett., 23, 3-4, 301-306 (1983).

Максимальная точность определена по ненулевым высшим производным X-i" = - Ю- - Ю_1 - ю2, x2" = —ю_2 - ro3(1-x!-x2)2 + 2ra3x2(1-x1-x2), x| ' ' = 0, x2"'

© В. Х. Федотов - канд. хим. наук, доц. каф. информационных систем ЧувГУ, [email protected]; Н. И.Кольцов -проф. каф. физической химии и ВМС ЧувГУ, [email protected].

д-р хим. наук,