Контроль устойчивости метода Ческино второго порядка точности Текст научной статьи по специальности «Математика»

УДК 519.622

© А.Е. Новиков, Е.А. Новиков, Л.В. Кнауб КОНТРОЛЬ УСТОЙЧИВОСТИ МЕТОДА ЧЕСКИНО ВТОРОГО ПОРЯДКА ТОЧНОСТИ1

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

Ключевые слова: явные методы, контроль точности и устойчивости, жесткие задачи.

© A.E. Novikov, E.A. Novikov, L.V. Knaub STABILITY CONTROL FOR A SECOND ORDER CESCHINO'S METHOD

The inequation for stability control for a Ceschino’s method of second order is constructed. Calculation results of stiff problems modeling are given, they confirm a growth of algorithm efficiency at the expense of stability control.

Keywords: explicit methods, control of accuracy and stability, stiff problems.

Введение

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

Обычно алгоритм управления шагом интегрирования строится на контроле точности численной схемы. Это естественно, потому что основным критерием является точность нахождения решения. Однако при применении алгоритмов интегрирования на основе явных формул для решения жестких задач данный подход приводит к потере эффективности и надежности [4]. Это связано с тем, что на участке установления противоречие между точностью и устойчивостью приводит к большому количеству повторных вычислений решения, а шаг выбирается значительно меньше максимально допустимого. Этого можно избежать дополнительным контролем устойчивости численной схемы. В настоящее время можно выделить два подхода к контролю устойчивости [3, 5]. Первый связан с оценкой максимального собственного числа матрицы Якоби через ее норму с последующим контролем (наряду с точностью) неравенства h\\fy\\<D [5], где положительная постоянная D связана с размером области устойчивости. Ясно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числа /-ГТ1;|Х матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства /7||/.|Т1;|Х||</.) [3]. Во всех рассмотренных ситуациях такая оценка не приводит к увеличению вычислительных затрат [4].

Здесь с применением предложенного в [6] способа оценки максимального собственного числа матрицы Якоби построено неравенство для контроля устойчивости метода Ческино второго порядка точности [7]. Приведены результаты расчетов, подтверждающие повышение эффективности алгоритма интегрирования за счет контроля устойчивости.

1. Метод Ческино

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений

У = fit-у), Ж)=y^h,^<tk, (1)

1Работа поддержана РФФИ (проекты 11-01-00106 и 11-01-00224)

138

где у и/- ^-мерные вещественные вектор-функции, ^ - независимая переменная, которая изменяется на заданном интервале [^0, 4]. Для решения (1) будем использовать явные формулы типа Рунге-Кутта следующего вида

Уп+1 =Уп+^Рт,к,

К = ¥ (*п,Уп) ,К=¥ (С + ^Куп + ^А), (2)

К = ¥ (К + \к>у» + \к2) -къ = ¥ <Х + к,уп+кг- 2 к2 +2къ),

где к - шаг интегрирования, кь 1</<4, - стадии метода, ртЬ 1</<4, - числовые коэффициенты, т -

порядок точности метода. При значениях коэффициентов

Ри ~ ^ ’ Р22 ~ —^, р2Ъ ~ 2, Р24 — 0 (3)

схема (2), (3) имеет второй порядок точности [7].

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

Согласно [7], численная формула (2) с коэффициентами

1 2 1

=7’ Р*2 =0,Аз =-,Рн =~ (4)

6 3 6

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

оценку ошибки 8И вида

4

3п=^Р^~Р2,)к,- (5)

1=1

В результате для контроля точности вычислений применяется неравенство

(6)

где ||-|| - некоторая норма в IV. 8 - требуемая точность расчетов. Имеет место соотношение дп=0(к3), шаг кас по точности выбирается по формуле hac=qh, где значение q находится из уравнения д3||8„||=8. Если д<1, то происходит повторное вычисление решения (возврат) с шагом к, равным qk. В противном случае вычисляется приближенное решение, а прогнозируемый шаг вычисляется по формуле

k„+1=qk. Неравенство (6) хорошо зарекомендовало себя при решении многих практических задач, и

ниже будет использоваться здесь. В дальнейшем алгоритм переменного шага на основе схемы (2) с коэффициентами (3), (4) и неравенством для контроля точности (6) назовем СЕ8СН42.

3. Контроль устойчивости численной схемы

Построим неравенство для контроля устойчивости схемы (2). Для этого применим (2) для решения линейной задачи у'=Ау с постоянной матрицей А. Первые три стадии к\, к2 и £3 схемы (2), применительно к данной задаче, имеют вид

к, =Хуп ,к2=(Х + ±Х2)уп ,к3=(Х + ±Х2 + ,

где Х=кА. Нетрудно видеть, что имеют место соотношения

1 , 1 ,

к\ ~ 2к2 +к3= —Хуп, 0.5 (к2-к1) = -х2уп.

О О

Теперь оценку максимального собственного числа матрицы Якоби системы (1) можно вычислить степенным методом. Введем обозначение

уя=2-тах|1№~2^2+^з)-||. (7)

п \{к2-кх),\ ] "

Тогда, согласно [6], для контроля устойчивости метода Ческино можно применять неравенство

где постоянное число В ограничивает интервал устойчивости.

Устойчивость методов типа Рунге-Кутта обычно исследуется на скалярном тестовом уравнении

139

г=1

/ = 4у, (9)

где X есть произвольное комплексное число, Яе(/.)<(). Смысл X - некоторое собственное число матрицы Якоби задачи (1). Так как X есть произвольное комплексное число с отрицательной вещественной частью, то для устойчивых задач все собственные числа матрицы Якоби системы (1) попадают в

(9). Применяя (2) с коэффициентами (3) для решения (9), получим, что функция устойчивости Q2(x) метода второго порядка точности имеет вид

1 о 1 ,

Оп (х) = 1 + х + — х' + — х~, х-Ьк.

Область устойчивости метода второго порядка приведена на рис. 1 (приведены линии уровня |<92(х)|=<7 при q равном 0.3, 07 и 1.0).

Рис. 1. Область устойчивости метода второго порядка точности

Из рис. 1 следует, что длина интервала устойчивости метода второго порядка точности (2) приблизительно равна 2. Функция устойчивости Q4(x) метода четвертого порядка точности имеет вид

ОЛх) = 1 + х + — х1 + — х3 + — х4, х = Ы.

~4 2 6 24

Область устойчивости метода четвертого порядка приведена на рис. 2. Интервал устойчивости метода четвертого порядка приблизительно равен 2.8. Поэтому в неравенстве (8) положим 0=2. Учитывая, что у„=0(И), шаг к по устойчивости можно выбирать по формуле к*=гк, где г вычисляется из равенства гу„=1).

Рис. 2. Область устойчивости метода четвертого порядка точности

Оценка (7) является грубой, потому что вовсе не обязательно максимальное собственное число сильно отделено от остальных, в степенном методе применяется мало итераций, и дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг вычисляется по

формуле

йя+1 =тах[йя,тт(/20С,/О], (10)

где кп есть последний успешный шаг интегрирования. Отметим, что формула (10) применяется для прогноза величины шага интегрирования кп+1 после успешного вычисления решения c предыдущим шагом кп, и поэтому фактически не приводит к увеличению вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого

140

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

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

4. Дифференциальные уравнения химической кинетики

Кинетическая схема химической реакции состоит из элементарных стадий, записываемых в виде

ыг ыг

1=1 1=1

где xг■, 1</<Ыг - реагенты; Nг и N - соответственно число реагентов и число стадий в реакции; а// и Ру, 1</<Ыг, 1</<Ы - стехиометрические коэффициенты. Для каждой элементарной реакции заданы константы скоростей стадии к/, 1</<Ы Процессу (11) в рамках сосредоточенной модели изотермического реактора постоянного объема соответствует система обыкновенных дифференциальных уравнений

С' = АТУ (12)

с заданным начальным условием С(0)=С0. Здесь Ат - стехиометрическая матрица, С и V - соответственно вектор концентраций реагентов и скоростей стадий. В случае протекания реакции в неизотермических условиях к данной системе добавляется уравнение теплового баланса

Г =

Є'Г-а(Г-Г01)

с;с

где Т - температура смеси в реакторе, Т01 - температура стенок реактора, $ - вектор удельных теп-лот стадий, Су - вектор теплоемкостей реагентов, а=ш,/г, а - коэффициент теплоотдачи, 5 и г - площадь поверхности и объем реактора соответственно. Верхний индекс Т у векторов $ и Су означает транспонирование. Теплоемкости реагентов и коэффициент теплоотдачи могут быть функциями концентраций реагентов сі, 1</<^г, а а также может еще зависеть от температуры.

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

, г 1

С = А V + —(С -С),

© '

где Ср - вектор концентраций реагентов на входе в реактор, 0 - время пребывания смеси в реакторе, 0=г/и, и - объемная скорость течения смеси через реактор. При протекании реакции в неизотермических условиях эта система дополняется уравнением теплового баланса

<?У-а{Т-Тт) 1

СТуС 0 02

где Т02 - температура смеси на входе в реактор. Температура реагирующей смеси может задаваться в виде функции времени ґ и концентраций сі, 1<і<^г, т.е. Т=Т(ґ,С).

5. Численное моделирование пиролиза этана

Пиролиз этана в отсутствии кислорода описывается небольшой последовательностью стадий. Механизм пиролиза этана неоднократно обсуждался в литературе. Здесь принята схема реакции, предложенная и исследованная в работе [8]

С2Я6 -> СЯ3+СЯ3,

СН3 + С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 = 3.69-103, к4 = 3.66-105, £5 = 1.62-107.

Обозначим концентрации реагентов следующим образом:

С1 = [С2Н6], с2 = [СН3], с3 = [СН4], с4 = [С2Н5], с5 = [С2Н4], с6 = [Н], с7 = [Н2], с8 = [С4Н10].

Тогда соответствующая система состоит из 8 обыкновенных дифференциальных уравнений вида

с[ — —к^с^ — к2схс2 — к4схс6, с' = 2^С[ — к2схс2, с’3 — к2схс2,

С4 к2схс2 к3сА -\-к4с^с^ 2к5с4. (13)

с' = к3с4, с' = к3с4 - к4схс6, с7 = к4схс6, с8 = к5с4 .

Начальная концентрация этана с1=[С2И6] равна 0.14, для остальных реагентов концентрации равны нулю.

Расчеты осуществлялись на РС Ше1^) Соге(ТМ) i7-3770S CPU@3.10GHz с двойной точностью с задаваемой точностью е=10-2. Эффективность алгоритмов интегрирования оценивалась по числу вычислений правой части if задачи (13) на интервале интегрирования. Численное решение осуществлялось на промежутке [0,0.26] с начальным шагом Ь=10-5. Данная задача удовлетворяет "классическому” определению жесткости. В начале интервала интегрирования наблюдается переходный участок (сотые доли секунды), а затем происходит медленное установление.

Сравнение эффективности построенного алгоритма (СЕ8СИ4281;) проводилось с известным методом Мерсона (МЕЯ80К) [9]. Для обоих методов фактическая точность не хуже задаваемой точности. Алгоритму СЕ8СИ4281 для нахождения решения потребовалось 20 403 вычислений правой части задачи (13), а для метода Мерсона if=25 796.

Заключение

Из результатов расчетов можно сделать следующие выводы. Во-первых, построенный алгоритм интегрирования второго порядка с контролем точности вычислений и устойчивости численной схемы можно применять для решения достаточно жестких задач. Во-вторых, по вычислительным затратам алгоритм СЕ8СИ4281 эффективнее метода Мерсона. Это является следствием контроля устойчивости численной схемы. Представляется, что при достаточно большой размерности задачи (1) метод СЕ8СИ4281 может конкурировать с неявными методами на задачах умеренной жесткости, потому что в нем не обращается матрица Якоби. При решении двенадцати тестовых задач [2] и десяти примеров

[10] преимущество алгоритма СЕ8СИ4281 значительно выше. В частности, при моделировании модифицированного орегонатора [11] преимущество СЕ8СИ4281 по сравнению с методом Мерсона более чем в полтора раза.

Литература

1. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений.

Нежесткие задачи. - М.: Мир, 1990. - 512 с.

2. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и

дифференциально-алгебраические задачи. - М.: Мир, 1999. - 685 с.

3. Новиков Е.А. Явные методы для жестких систем. - Новосибирск: Наука, 1997. - 197 с.

4. Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жестких гибридных систем. -

Новосибирск: Изд-во НГТУ. 2012. - 450 с.

5. Shampine L.M. Implementation of Rosenbrock methods // ACM Transaction on Mathematical Software. - 1982. №5. P. 93-113.

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

7. Ceschino F., Kuntzman J. Numerical solution of initial value problems. - New Jersey: Prentice-Hall, Englewood Clis, 1966. - 318 p.

8. Kulich D.M., Taylor J.E. Mathematical simulation of the oxygen ethane reaction // J. Chem. Kinetic. 1975. V. 8. P. 89-97.

9. Merson R.H. An operational methods for integration processes // Australia: Proc. of Symp. on Data Processing. 1957. P. 329-331.

10. Enright W.H., Hull T.E. Comparing numerical methods for the solutions of systems of ODE’s // BIT. - 1975. - №15. - P. 10-48.

11.Новиков Е.А. Численное моделирование модифицированного орегонатора (2,1)-методом решения жестких задач // Вычислительные методы и программирование. 2010. -Т.11. - №2. - С. 123-130.

Новиков Антон Евгеньевич, преподаватель кафедры Сибирского федерального университета, тел. +7(391)2495979. E-mail: aenovikov@bk .ru.

Новиков Евгений Александрович, главный научный сотрудник Института вычислительного моделирования СО РАН, тел. +7(391)2494724. E-mail: novikov@icm .krasn.ru.

Кнауб Людмила Владимировна, доцент кафедры Сибирского федерального университета, тел. +7(391)2912213. E-mail: aenovikov@bk.ru

Novikov Anton Evgenevich, lecturer, Siberian Federal University, ph. +7(391)2495979.

E-mail: aenovikov@bk .ru.

Novikov Evgeny Alexandrovich, chief researcher, Institute of Calculation Modeling SB RAS, ph. +7(391)2494724. E-mail: novikov@icm .krasn.ru.

Knaub Lyudmila Vladimirovna, associate professor, Siberian Federal University, ph. +7(391)2912213. E-mail: aenovikov@bk.ru