CC BY
45
УДК 519.622
Е. А. Новиков
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПИРОЛИЗА ЭТАНА ЯВНЫМ МЕТОДОМ ТРЕТЬЕГО ПОРЯДКА ТОЧНОСТИ1
Аннотация. Получены коэффициенты явного трехстадийного метода типа Рунге - Кутта. Построены неравенства для контроля точности вычислений и устойчивости численной схемы. Результаты моделирования пиролиза этана демонстрируют повышение эффективности за счет дополнительного контроля устойчивости.
Ключевые слова: жесткая задача, явный метод, контроль точности и устойчивости, пиролиз этана.
Abstract. Coefficients of explicit three-stage Runge - Kutta method have been obtained. The inequalities for exactness of calculations control and stability control of numerical scheme have been developed. Result numerical modeling of ethane pyrolysis demonstrate an efficiency increase with an additional stability control.
Keywords: stiff problem, explicit method, control accuracy and stability, pyrolysis of ethane.
Введение
Во многих приложениях возникает необходимость решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Обычно для численного решения таких задач применяют алгоритмы на основе неявных или полуявных численных формул вследствие их хороших свойств устойчивости. В данных методах используется декомпозиция матрицы Якоби с выбором главного элемента по строке или столбцу, а иногда и по всей матрице. При большой размерности исходной системы это отдельная трудоемкая задача. Если элементы матрицы Якоби носят нерегулярный характер, то получение данной матрицы и составление подпрограммы ее нахождения требуют от вычислителя больших затрат времени. Это характерно, например, для дифференциальных уравнений химической кинетики. При численном определении данной матрицы возникает проблема с выбором шага численного дифференцирования. В такой ситуации предпочтительнее применять алгоритмы на основе явных численных формул, если жесткость задачи позволяет за разумное время получить приближение к решению [1].
Современные алгоритмы на основе явных методов в большинстве своем не приспособлены для решения жестких задач по следующей причине. Обычно алгоритм управления шагом интегрирования строится на контроле точности численной схемы. Это естественно - основным критерием является точность нахождения решения. Однако при применении таких алгоритмов для решения жестких задач этот подход приводит к потере эффективности и надежности, потому что на участке установления вследствие противоречивости требований точности и устойчивости шаг интегрирования раскачивается. В лучшем случае это приводит к большому количеству повторных вычислений решения, а шаг выбирается значительно меньше допустимого. Этого
1 Работа поддержана грантами РФФИ № 08-01-00621 и Президента НШ-3431.2008.9
можно избежать, если наряду с точностью контролировать устойчивость численной схемы.
В настоящее время можно выделить два подхода к контролю устойчивости [2-3]. Первый способ связан с оценкой максимального собственного числа матрицы Якоби Уу через ее норму с последующим контролем (наряду с контролем точности) неравенства ЩУУ\\ < D [2], где h есть шаг интегрирования, а положительная постоянная D зависит от размера области устойчивости метода. Ясно, что для явных методов, в которых матрица ЯкобиУу не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к значительному увеличению вычислительных затрат.
Второй подход основан на оценке максимального собственного числа ^тах матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства ^^тах\ < D [3]. Во всех рассмотренных ситуациях такая оценка фактически не приводит к увеличению вычислительных затрат [1, 3].
В данной работе построен алгоритм интегрирования переменного шага на основе трехстадийной схемы типа Рунге - Кутта третьего порядка с контролем точности вычислений и устойчивости численной схемы. На примере моделирования пиролиза этана продемонстрировано повышение эффективности расчетов за счет дополнительного контроля устойчивости.
1. Численная схема
Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений
У = У(У), у(^о) = Уo, ^ ^ t ^ ^ , (1)
рассмотрим явный трехстадийный метод типа Рунге - Кутта вида
Уп+1 = Уп + Р1Ч1 + Р2Ч2 + Р3Ч3,
Ч = ¥ (Уп ), Ь = ¥ (Уп + Р21Ч1),
Ч3 = ¥(уп + Рз1Ч1 + Рз2 Ч2) , (2)
где у и У - вещественные Л-мерные вектор-функции; t - независимая переменная; h - шаг интегрирования; Чь Ч2 и Ч3 - стадии метода; р1, р2, р3, р2Ь р31 и р32 - числовые коэффициенты, определяющие свойства точности и устойчивости (2). В случае неавтономной задачи
У' = у^,У), У(t0) = Уо , t0 ^t ^ Ь ,
схема (2) записывается в виде
Уп+1 = Уп + Р1Ч1 + Р2Ч2 + Р3Ч3 ,
Ч1 = ¥ ^п>Уп ),
Ч2 = ¥ (tn + Р21^ уп + Р21Ч1) ,
Ч3 = ¥(tn +[Р31 + Р32]^ уп +Р31Ч1 +Р32Ч2) .
Ниже для сокращения выкладок будем рассматривать (1). Однако построенные далее методы можно применять для решения неавтономных задач.
Получим соотношения на коэффициенты метода (2) третьего порядка точности. Для этого разложим стадии Ч1, Ч2 и Ч3 в ряды Тейлора по степеням h до членов с '4 включительно и подставим в первую формулу (2). В результате получим
уп+1 = уп + (Р1 + Р2 + p3)h/n + [Р21Р2 + (Р31 + Р31) Р3 ]h 2упуп +
+
Р21Р32Р3уп уп + 2[Р21Р2 + (Р31 + Р32)2Р3
12
Р21Р32 Р3 Гп^/п + Р21(Р31 + Р32)Р32 Р3уппупуп +
+ 6(1 ^2 + (Р31 + Р32)3Р3
+ O(h5),
(3)
где элементарные дифференциалы вычислены на приближенном решении уп,
т.е. уп = у(Уп), =д/ (Уп)/ЭУ, = д2у(Уп)/ЭУ2 и =д3у(Уп)/Эу3. Точ-
ное решение у(^) в окрестности точки tn имеет вид
у (tn+l) = у «п) + hf + 0,5h2 Г/ +1 h3 Г У' 2У + УУ2
24
У '3У + УГУ2 + 3УУУ2 + УУ3
+ 0('5);
(4)
где элементарные дифференциалы вычислены на точном решении у(^).
Сравнивая полученные ряды для приближенного (3) и точного (4) решений до членов с ' включительно при условии уп = у(^), запишем условия третьего порядка точности схемы (2):
Р1 + Р2 + Р3 = 1,
р21Р2 + (р31 + Р32) Р3 = 0,5 ,
р21Р2 + (Р31 + Р32)2Р3 = 1/3 ,
Р21Р32 Р3 =1/6
(5)
В предположении уп = y(tn) локальную ошибку 5п+1 схемы (2) можно вычислить по формуле 5п+1 = у(^+1) - уп+1. Учитывая представления уп+1 и y(tn+1) в виде рядов Тейлора (3) и (4), запишем
п+1
24
'А\^У,ЪУ+1 Р21^32Р3 IУУУ2 +
24 2
--Р21(Р31 +Р32)Р32 р3
X
X/ff2 +
ТГ-6Р21Р2 -7(Р31 +Р32)3Р3 24 6 6
ГУЪ\ + о('5). (6)
Далее получим коэффициенты численной схемы (2) третьего порядка точности.
2. Исследование условий порядка
В нелинейной системе алгебраических уравнений (5) два свободных коэффициента. Исследуем три варианта.
Вариант 1. Положим р21 = Рзі + Рзг и р31 = р32. Это означает, что приращения к2 и к3 будут вычислены в одной и той же точке їп + р21й, причем вклад к1 и к2 при определении к3 учитывается одинаково. Тогда нелинейную систему (5) можно переписать в виде
Из второго и третьего уравнений данной системы имеем р21 = 2/3. Из соотношений р21 = р31 + р32 и р31 = р32 запишем р31 = р32 = 1/3. Из четвертого уравнения системы получим р3 = 3/4. Из равенства р2 + р3 = 3/4 имеем р2 = 0. Наконец, из первого соотношения системы получим р1 = 1/4. В результате коэффициенты схемы (2) определяются однозначно и имеют вид
При данных соотношениях локальную ошибку 5п+1 схемы (2) можно записать следующим образом:
Вариант 2. Минимизируем локальную ошибку (6). Для этого, учитывая вид (6), вместо (5) рассмотрим следующую расширенную нелинейную систему алгебраических уравнений:
1) р1 + р2 + р3 =1;
2) Р21(р2 + р3) = -2
3) р21(р2 + р3) = -
4) Р21Р32 р3 = 6 .
6
(7)
8и+1 = -^к4 \9/3/ + //3 - 3/Я2 + 3///2 ] + 0(Н5).
216
1) р1 + р2 + р3 =1;
2) Р21р2 + (Р31 + Р32)р3 = 2 ;
3) $21р2 + (Р31 +Р32)2 р3 = 3
4) Р21Р32 р3 = Т ;
6
5) Р21Р32р3 = 12 ;
6) Р21(Р31 +Р32)Р32 р3 = 3.
Исследуем совместность данной системы. При 3р21/2 = р31 + р32 два последних уравнения совпадают. Из четвертого и пятого соотношений имеем р21 = 1/2. Из второго и третьего равенств получим р2 = 1/3 и р3 = 4/9. Из первого уравнения запишем р1 = 2/9, а из четвертого имеем р32 = 3/4. Наконец, из соотношения р31 + р32 = 3/4 запишем р31 = 0. В результате коэффициенты метода (2) с минимальной локальной ошибкой можно записать в виде
1 3 2 14
р21 = 2, Р31 = 0 , Р32 = 4, р1 = 9, р2 = 3 , р3 = 9 . (8)
При данных соотношениях локальную ошибку 5п+1 схемы (2) можно представить следующим образом:
1
5
п+1 ^
п+1 288
12 У ,ЪУ - УУ3
+ О('5).
При использовании (2) с наборами коэффициентов (7) или (8) ни одна стадия не вычисляется в точке ^п+1. При быстром изменении решения это может приводить к понижению эффективности расчетов.
Вариант 3. Положим р21 = 1/2 и р31 + р32 = 1. Тогда на каждом шаге приращения к1, к2 и к3 вычисляются соответственно в точках ^п, ^п + 0,5' и Iп + '. В этом случае условия третьего порядка записываются в виде
1) Р1 + Р2 + Р3 =1;
1 1
2) 2 р2 + р3 =
04 1 1
3) - Р2 + Р3 = 3;
4) Р32 р3 = 3 .
Из второго и третьего равенств данной системы запишем Р2 = 2/3 и р3 = 1/6. Из первого и последнего уравнений имеем р1 = 1/6 и р32 = 2. Из равенства р31 + р32 = 1 следует р31 = -1. В результате коэффициенты метода (2) имеют вид
1 12 1
Р21 = ^ Р31 =-1, Р32 = ^ р1 = ^ р2 = -, р3 = -. (9)
2 6 3 6
При данных соотношениях локальную ошибку 5п+1 схемы (2) можно записать следующим образом:
5п+1 =
1 '4
72
3УвУ - зfff2 - УУ3 + О(¥).
5
3. Контроль точности вычислений
Построим неравенство для контроля точности вычислений метода третьего порядка. Для этого рассмотрим вспомогательную схему
Уп+1Д = Уп + г1к1 + г2к2 ,
где к1 и к2 определены в (2). Потребуем, чтобы данный метод имел второй порядок точности. Разложение приближенного решения уп+11 в виде ряда Тейлора по степеням ' имеет вид
уп+1,1 = уп + (г1 + г2)Г + Р21г2'2/п/п + О('3) .
Сравнивая ряды Тейлора для точного у^О и приближенного уп+11 решений, видим, что требование второго порядка точности будет выполнено, если
П + г2 = 1, Р21г2 = 0,5.
Отсюда получим
г2 = 0,5/Р21, п = 1 - Г2,
где значение р21 определено в (7), (8) или (9). Теперь с помощью идеи вложенных методов оценку аналога глобальной ошибки еп,3 метода третьего порядка точности можно вычислить по формуле [1]
еп,3 = Уп+1 - Уп+1,1 =(Р1 - г1)к1 +(Р2 - п2)к2 + Рзк3 .
Тогда неравенство для контроля точности вычислений имеет вид
11(Р1 - п1)к1 + (Р2 - г2)к2 + Р3к3 У< е ,
где ||^|| - некоторая норма в В1; е - требуемая точность интегрирования. В конкретных расчетах применялся метод (2) с коэффициентами (9) как более надежный. Тогда неравенство для контроля точности имеет вид
|| к1 - 2к2 + к3||< 6е. (10)
3. Контроль устойчивости численной схемы
Теперь построим неравенство для контроля устойчивости численной формулы (2) предложенным в [1] способом. Запишем стадии к1, к2 и к3 применительно к задаче у' = Ау, где А - матрица с постоянными коэффициентами. В результате получим
к1 = Хуп , к2 = (Х + Р21Х2)уп ,
к3 =
Х + (Р31 + Р32)Х 2 + Р21Р32Х 3
уп ,
где X = 'А.
Найдем коэффициенты й1, й2 и й3 из условия
3
+ ^2к2 + <^зкз = X Уп .
Данное требование будет выполнено, если
й1 = (Р31 + Р32 - Р21)/ й , й2 = _(Р31 + Р32)/ й , й3 = Р21 / й , где й = р221р32. Нетрудно видеть также, что
Р21 (к2 - к1) = Х 2 уп .
Тогда согласно [1] оценку максимального собственного числа vn,3 = h|Amax| матрицы Якоби системы (1) можно вычислить по формуле
vn 3 = P21 max (| d( + d2k2 + d^k1^ \ l\ k2 -k{
1<i< N'
Интервал устойчивости численной схемы (2) приблизительно равен 2,5. Поэтому для ее контроля устойчивости можно применять неравенство Vn,3 < 2,5. Полученная оценка является грубой, потому что:
1) вовсе не обязательно максимальное собственное число сильно отделено от остальных;
2) в степенном методе применяется мало итераций;
3) дополнительные искажения вносит нелинейность задачи (1).
Поэтому здесь контроль устойчивости используется как ограничитель
на размер шага интегрирования. В результате прогнозируемый шаг hn+1 будем вычислять следующим образом.
Новый шаг hac по точности определим по формуле hac = q1hn, где hn есть последний успешный шаг интегрирования, а q1, учитывая соотношение sn,3 = O(h3n), задается уравнением q31||en,3|| = s. Шаг hs по устойчивости зададим формулой hst = q2hn, где q2, учитывая равенство vn 3 = O(hn), определяется из уравнения q2vn3 = 2,5. В результате шаг hn+1 вычисляется по формуле
hn+1 = max hn, min(hac, hst)
Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг hn. В результате для выбора шага и предлагается данная формула. Она позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Собственно говоря, именно наличие этого участка существенно ограничивает возможности применения явных методов для решения жестких задач.
4. Численное моделирование пиролиз этана
Расчеты проводились на Intel(R) Core 2 Quad CPU с двойной точностью. В конкретных расчетах левая часть неравенства для контроля точности (10) вычислялась по формуле
|| k1 - 2k2 + кз || = max \ k1 - 2k2 + k3 \ /(\ yln \ +r)
1<i<N L
где i - номер компоненты; r - положительный параметр.
Если по i-й компоненте решения выполняется неравенство |yn!| < r, то контролируется абсолютная ошибка re, в противном случае - относительная ошибка s. В расчетах параметр r выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой.
Пиролиз этана в отсутствие кислорода описывается небольшой последовательностью стадий. Механизм пиролиза этана неоднократно обсу-
ждался в литературе. Здесь принята схема реакции, предложенная и исследованная в [4]
С2Н6 — СНз + СНз,
СНз + С2Н6 —— С2Н4 + С2Н5,
С2Н5 — С2Н4 + Н,
Н + С2Н6 — Н2 + С2Н5,
С2Н5 + С2Н5 — С4Н10, где константы скоростей стадий имеют вид
к = 1,34 ■ 10-5, к2 = 3,73 ■ 102, кз = 3,69 ■ 103, к4 = 3,66 ■ 105, к5 = 1,62 ■ 107.
Обозначим концентрации реагентов следующим образом:
С = [С2Н6], С2 = [СНз], Сз = [СН4], С4 = [С2Н5],
С5 = [С2Н4], С6 = [Н], С 7 = [Н2], С8 = [С4Н10].
Соответствующую систему дифференциальных уравнений можно получить с применением алгоритма, описанного в [5]. Данная система состоит из восьми уравнений и имеет вид
С1 = _к1С1 -к2С1С2 -к4С1С6 , С2 = 2к1С1 - к2С1С2 , С3 = к2С1С2 ,
С4 = к2С1С2 - кзС4 + к4С1С6 - 2к5С2 , С5 = кзС4,
С6 = к3С4 - к4С1С6 , с7 = к4С1С6 , с8 = к5С4 . (11)
Начальная концентрация этана С1 = [С2Н6] равна 0,14, для остальных реагентов начальные концентрации равны нулю.
Расчеты осуществлялись с точностью е = 10-4, при которой наиболее эффективны методы третьего порядка. Эффективность алгоритмов интегрирования оценивалась по числу вычислений правой части гУ задачи (11) на интервале интегрирования. Численное решение осуществлялось на промежутке [0; 0,26] с начальным шагом ' = 10-5. Данная задача удовлетворяет «классическому» определению жесткости. В начале интервала интегрирования наблюдается переходный участок (сотые доли секунды), а затем происходит медленное установление.
Сравнение эффективности алгоритма интегрирования без контроля устойчивости (ИК3) и с контролем устойчивости (ИК38Т) проводилось известным методом Мерсона [6]. Для всех методов фактическая точность не хуже задаваемой точности. Алгоритму КК3 для нахождения решения потребовалось 19 790 вычислений правой части задачи (11), для алгоритма ИК38Т гУ = 17 004, а для метода Мерсона 1У = 26 876. Таким образом, на умеренно жестких задачах построенный алгоритм с контролем точности и устойчивости примерно в полтора раза эффективнее метода Мерсона.
Заключение
Из результатов расчетов можно сделать следующие выводы. Во-первых, построенный алгоритм интегрирования третьего порядка с контро-
лем точности вычислений и устойчивости численной схемы можно применять для решения достаточно жестких задач. Во-вторых, по вычислительным затратам алгоритм RK3ST эффективнее метода Мерсона примерно в 1,5 раза. Это является следствием контроля устойчивости численной схемы. Представляется, что при достаточно большой размерности задачи (11) метод RK3ST может конкурировать с неявными методами на задачах умеренной жесткости, потому что в нем не обращается матрица Якоби.
Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f . Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки.
Список литературы
1. Новиков, Е. А. Явные методы для жестких систем / Е. А. Новиков. - Новосибирск : Наука, 1997.
2. Shampine, L. M. Implementation of Rosenbrock methods / L. M. Shampine // ACM Transaction on Mathematical Software. - 1982. - V. 8. - № 5. - P. 93-113.
3. Новиков, В. А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных дифференциальных уравнений /В. А. Новиков, Е. А. Новиков // ДАН СССР. - 1984. - Т. 277. - № 5. - С. 1058-1062.
4. Kulich, D. M. Mathematical simulation of the oxygen ethane reaction / D. M. Ku-lich, J.E. Taylor // J. Chem. Kinet. - 1975. - V. 8. - P. 89-97.
5. Новиков, Е. А. Комплекс программ моделирования кинетики сложных реакций / Е. А. Новиков, Ю. А. Шитов, В. И. Бабушок, Д. В. Марьин // Прямые и обратные задачи в химической кинетике. - Новосибирск : Наука, 1993. - С. 22-38.
6. Merson, R. H. An operational methods for integration processes / R. H. Merson // Proc. of Symposium on Data Processing. Weapons Research Establishment. - Australia : Salisbury, - 1957. - P. 329-331.
Новиков Евгений Александрович
доктор физико-математических наук, профессор, главный научный сотрудник, Институт вычислительного моделирования СО РАН (Красноярск)
E-mail: novikov@icm.krasn.ru
Novikov Evgeny Alexandrovich Doctor of Physical and Mathematical Sciences, chief researcher, Institute of computational modeling of the Russian Academy of Sciences (Krasnoyarsk)
УДК 519.622 Новиков, Е. А.
Численное моделирование пиролиза этана явным методом третьего порядка точности / Е. А. Новиков // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 4 (16). -С. 64-72.
