Алгоритм переменного порядка на основе стадий метода Ческино

Тип работы:
Реферат
Предмет:
Физико-математические науки


Узнать стоимость

Детальная информация о работе

Выдержка из работы

© Е.А. НОВИКОВ, А.А. ЗАХАРОВ
novikov@icm. krasn. ru, azaharov@utmn. ru
УДК 519. 622
АЛГОРИТМ ПЕРЕМЕННОГО ПОРЯДКА НА ОСНОВЕ СТАДИЙ МЕТОДА ЧЕСКИНО*
АННОТАЦИЯ. В работе исследуются методы численного решения жестких задач большой размерности. С помощью оценки максимального собственного числа матрицы Якоби построено неравенство для контроля устойчивости численной схемы Ческино второго порядка точности. Для интегрирования с переменным шагом предложена формула, позволяющая прогнозировать следующий шаг по времени. На основе этой формулы разработан метод первого порядка точности с расширенной областью устойчивости. Этот метод позволяет стабилизировать поведение шага интегрирования на участке установления решения, где именно устойчивость играет определяющую роль. Тем самым удается снять ограничения на возможность применения явных методов для решения жестких задач. Сформулирован алгоритм численного решения жестких задач переменного порядка, использующий неравномерный шаг по времени с дополнительным контролем устойчивости численной схемы интегрирования. Приведены результаты расчетов жестких задач, связанных с численным моделированием пиролиза этана. Получено подтверждение повышения эффективности расчетов за счет переменного порядка.
SUMMARY. This paper investigates methods for the numerical solution of stiff problems with large dimension. Using the estimation of the largest eigenvalue of the Jacobi matrix, we create an inequality in order to control the stability of a Cescino numerical scheme with second-order accuracy. In order to integrate a variable step, we propose a formula which allows one to predict the next step in time. On the basis of this formula, we develop a method with first-order accuracy with extended stability range. This method allows stabilized behaviour of step integration at the stage of solution exactly where stability plays a crucial role. This makes it possible to remove restrictions on the possibility of using explicit methods for solving stiff problems. We formulate an algorithm for the numerical solution of stiff problems of variable order, which uses the irregular step in time with an additional control of stability of the numerical integration scheme. This paper demonstrates solutions of stiff problems associated with numerical simulations of ethane pyrolysis, which confirm an increase in efficiency due to the use of variable order.
КЛЮЧЕВЫЕ СЛОВА. Явные методы, контроль точности и устойчивости, жесткие задачи.
KEY WORDS. Explicit methods, accuracy and stability control, stiff problems.
* Работа выполнена при поддержке РФФИ (проект 11−01−106)
Введение. При численном решении жестких задач большой размерности возникает необходимость использования алгоритмов на основе явных методов [1−3]. Методы интегрирования на основе неявных или полуявных численных схем, как правило, используют обращение матрицы Якоби [2]. В данном случае это отдельная трудоемкая задача. В такой ситуации предпочтительнее применять алгоритмы на основе явных формул, если жесткость задачи позволяет за разумное время получить приближение к решению [3].
Алгоритм управления шагом интегрирования строится на контроле точности численной схемы. Основным критерием является точность нахождения решения. Однако при применении алгоритмов интегрирования на основе явных формул для решения жестких задач данный подход приводит к потере эффективности и надежности [4−5]. Это связано с тем, что на участке установления противоречие между точностью и устойчивостью приводит к большому количеству повторных вычислений решения, а шаг выбирается значительно меньше максимально допустимого. Этого можно избежать дополнительным контролем устойчивости численной схемы. Обычно используют два подхода к контролю устойчивости [5−7]. Первый связан с оценкой максимального собственного числа матрицы Якоби через ее норму с последующим контролем (наряду с точностью) неравенства h\fy\& lt-D [6], где положительная постоянная D связана с размером области устойчивости. Для явных методов, где матрица Якоби не участвует в вычислительном процессе. Это приводит к увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числашах матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства [7]. Во всех рассмотренных ситуациях такая оценка не
приводит к увеличению вычислительных затрат [3−5, 7].
Здесь с применением предложенного в [7] способа оценки максимального собственного числа матрицы Якоби построено неравенство для контроля устойчивости метода Ческино второго порядка точности [8]. Приведены результаты расчетов, подтверждающие повышение эффективности алгоритма интегрирования за счет контроля устойчивости.
Метод Ческино. Рассмотрим задачу Коши
где у и / --мерные вещественные вектор-функции, t — независимая переменная, которая изменяется на заданном интервале [^, tk]. Для решения (1) будем использовать явные формулы типа Рунге-Кутта
где h — шаг интегрирования,, 1& lt-г<-4, — стадии метода, рт1, 1& lt-г<-4, — числовые коэффициенты, т — порядок точности метода. При коэффициентах
У'-=ї(і, У), Ж)= Уо, іо& lt-і<- ік
(1)
yn*1_yn+pm1k1+pm2k1+pm3k3+pm4k4 ,
к^(іп, у"), к2=hf (іn+h/4, у"+к/4), кз=hf (іn+h/2, Уп+ к2/2), к4^/(іп+^ уп+ ^-2^+2^),
(2)
(3)
Контроль точности вычислений. Схема (2) с коэффициентами
р41=1/6, Р42=0, Р43=2/3, р44=1/6
имеет четвертый порядок. Тогда для контроля точности схемы второго порядка можно использовать оценку ошибки 5п2 вида
5″, 2=(Р41-Р21)к1+(Р42-Р22)к2+(Р43-Р23)к3+(Р44-Р24)к4-
В результате для контроля точности вычислений применяется неравенство \5″, 2\& lt-в, где \-\ - некоторая норма в Ям, в — требуемая точность расчетов. Учитывая 5п, 2=0(К3), шаг Кас по точности выбирается по формуле hac=qh, где значение q находится из уравнения q3\ 5п, 2\=в. Если q& lt-1, то происходит повторное вычисление решения (возврат) с шагом h, равным qh. В противном случае находится приближенное решение, а прогнозируемый шаг hn+^ вычисляется по формуле КnA=qК. Неравенство \5п, 2\& lt-в хорошо зарекомендовало себя при решении многих практических задач. Алгоритм переменного шага на основе схемы (2), (3) с неравенством для контроля точности \5п, 2\& lt-в назовем СЕБСН42.
Контроль устойчивости численной схемы. Построим неравенство для контроля устойчивости схемы (2). Для этого применим (2) для решения линейной задачи у'-=Ау с постоянной матрицей А. Первые три стадии к1, к2 и к3 применительно к данной задаче имеют вид
к=Хуп, к2=(Х+Х2/4)уп, к,=(Х+Х2/2+Х3/8)уп, где Х=КА. Нетрудно видеть, что имеют место соотношения к-2к2+к3=Х3уп/8, 0. 5(к2-к1)=Х2уп/8.
Теперь оценку максимального собственного числа матрицы Якоби системы
(1) можно вычислить степенным методом [3]. Введем обозначение
Уп=2-шах1& lt--<-Л,{ Кк-2к2+к3)г/((к2-к:)^) }. (4)
Тогда для контроля устойчивости метода Ческино можно использовать неравенство уп^, где число D ограничивает интервал устойчивости.
Устойчивость методов типа Рунге-Кутта обычно исследуется на скалярном тестовом уравнении у'-=ку, где есть произвольное комплексное число, Ке (^)& lt-0. Смысл — некоторое собственное число матрицы Якоби задачи (1). Применяя
(2), (3) для решения у'-=ку получим, что функция устойчивости Q2(x) метода второго порядка точности имеет вид
д2(х)=(1+Х+Х2/2+Х¾), х=Ь^, а функция устойчивости Q4(x) метода четвертого порядка:
Q4(x)=(1+ Х+Х2/2+Х3/6+Х4/24), х=Кк.
Интервал устойчивости метода второго порядка равен 2, а схемы четвертого порядка ~ 2.8. Поэтому в неравенстве уп^ положим D=2. Учитывая, что уп=0(К), шаг К^ по устойчивости можно выбирать по формуле К^=гК, где г вычисляется из равенства гуп=2.
Оценка (4) является грубой, потому как не обязательно, что максимальное собственное число сильно отделено от остальных собственных значений. В степенном методе применяется мало итераций и дополнительные искажения вносит нелинейность задачи (1). Поэтому контроль устойчивости используется
как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг вычисляется по формуле
Кп, 1=шах{Кп, тт[Кас, К5*]}, (5)
где Кп есть последний успешный шаг интегрирования. Отметим, что формула
(5) применяется для прогноза величины шага интегрирования Кп+1 после успешного вычисления решения с предыдущим шагом Кп, и поэтому фактически не приводит к увеличению вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет увеличен, т.к. не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг Кп. В результате для выбора шага и предлагается формула (5). Данная формула позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Именно наличие данного участка существенно ограничивает возможности применения явных методов для решения жестких задач.
В дальнейшем алгоритм переменного шага с дополнительным контролем устойчивости численной схемы будем называть CESCH42st. Данный алгоритм основан на численной формуле невысокого (второго) порядка точности, и поэтому он ориентирован на решение нежестких задач с небольшой точностью расчетов (порядка 1% и ниже), а также задач умеренной жесткости. Из результатов расчетов алгоритмом СЕ5СН42з1 следует, что фактическая точность вычислений на участке установления значительно выше задаваемой, потому что на данном участке старые ошибки подавляются за счет контроля устойчивости, а новые ошибки невелики за счет малости производных решения. В такой ситуации эффективнее проводить расчеты по методу низкого порядка с более широкой областью устойчивости.
Метод первого порядка. На основе стадий численной схемы (2) построим метод первого порядка точности вида
Упт=Уп+р1А+Р12к1+Ртк3+Рик4 (6)
с более широкой областью устойчивости. Для этого применим (6) для решения тестового уравнения у'-=& quot-ку. Получим yn+1=Q1(x)yn, где х=Кк, а функция устойчивости Q1(x) имеет вид:
Требование первого порядка точности численной формулы (7) означает выполнение соотношения р11+р12+р13+р14=1. Остальные коэффициенты ри применим для расширения области устойчивости. Условие устойчивости схемы (6) имеет вид Q1(x) & lt-1. Для построения метода с максимальным интервалом устойчивости рассмотрим многочлен Чебышева вида Г4(г)=8г4−8г2+1. Известно, что полином Т,(7.) наименее отклоняется от нуля при ге[-1,1]. Проведем замену переменных z=1−2x/у, при этом отрезок [у, 0] отображается на промежуток [-1,1]. В результате многочлен Т4^) записывается в виде
Для T4(x) неравенство Т^)^ выполняется на максимальном интервале [у, 0], у=-32 [3]. Сравнивая соотношения (7) и (8) при у=-32, получим коэффициенты ри, 1& lt-г<-4, метода (2) первого порядка точности с максимальным интервалом устойчивости:
Область устойчивости метода первого порядка (6), (9) по вещественной оси в 16 раз шире области устойчивости численной схемы (2), (3). Вычислительные затраты в методах первого и второго порядков совпадают. Поэтому для задач, в которых шаг ограничен в основном по устойчивости, предполагаемое теоретическое повышение эффективности в 16 раз.
Контроль точности и устойчивости метода (6), (9). В неравенстве для контроля точности будем применять оценку локальной ошибки
Тогда для контроля точности численной формулы (6), (9) можно применять неравенство 5п1& lt-в. Далее, так как интервал устойчивости численной схемы
(6), (9) ограничен числом 32, то для ее контроля устойчивости можно применять неравенство уп& lt-32, где уп определяется по формуле (4).
Алгоритм переменного порядка. Методы первого порядка с расширенными областями устойчивости эффективны на участках установления, где шаг ограничен по устойчивости. При высокой точности расчетов на переходных участках эффективнее будут методы (2), (3). Повышения эффективности можно достигнуть за счет применения каждого метода на том участке, где он наиболее эффективен. В качестве критерия переключения с метода на метод можно использовать неравенство для контроля устойчивости. При расчетах по методу (2), (3) переход на численную схему (6), (9) осуществляется при нарушении неравенства уп& lt-2. При расчетах методом первого порядка обратный переход происходит в случае выполнения уп& lt-2. Вычисления методом первого порядка сопровождаются дополнительным (наряду с точностью) контролем неравенства уп& lt-32, а шаг выбирается по формуле типа (5).
Дифференциальные уравнения химической кинетики. Кинетическая схема химической реакции состоит из элементарных стадий вида
где xi, 1& lt-г^К — реагенты- NR и N5 — соответственно, число реагентов и число стадий в реакции- а-/ и Р-/, 1& lt-г<- NR, 1& lt-/^5 — стехиометрические коэффициенты. Для каждой элементарной реакции заданы константы скоростей стадии 1 & lt-/<-N5.
Процессу (10) в рамках сосредоточенной модели изотермического реактора постоянного объема соответствует система обыкновенных дифференциальных уравнений С'-=АТУ с заданным начальным условием С (0)=С0. Здесь А Т — стехиометрическая матрица, С и V — соответственно, вектор концентраций реагентов и скоростей стадий. В случае протекания реакции в неизотермических условиях к данной системе добавляется уравнение теплового баланса T'-={QтV-а (Т-Т01)}/ {С^С}, где Т — температура смеси в реакторе, Т01 — температура стенок реактора, Qт — вектор удельных теплот стадий, С^ - вектор теплоемкостей реагентов,
а=аз/г, а — коэффициент теплоотдачи, в и г — площадь поверхности и объем реактора, соответственно. Верхний индекс Т у векторов Qт и С^ означает транспонирование. Теплоемкости реагентов и коэффициент теплоотдачи могут быть функциями концентраций реагентов с-, 1& lt-г^Я, а а может еще зависеть от температуры.
Если реакция протекает в изотермическом реакторе постоянного объема с обменом вещества (открытая система, реактор идеального смешения), система дифференциальных уравнений записывается в виде C'-=ATV+(Cp-C)/©, где Ср — вектор концентраций реагентов на входе в реактор, (c) — время пребывания смеси в реакторе, (c)=г/и, и — объемная скорость течения смеси через реактор. При протекании реакции в неизотермических условиях, эта система дополняется уравнением теплового баланса Т'-={^^-а (Т-Т01)}/(С^, ТС) — (Т-Т02)/©, где Т02 — температура смеси на входе в реактор. Температура реагирующей смеси может задаваться в виде функции времени t и концентраций реагентов сI, 1& lt-1<-МГ, то есть Т=Т (^С).
Алгоритм формирования уравнений химической кинетики [9]. Если стадия обратима, то скоростью стадии принимается разность скоростей прямого Ш+в и обратного W-s процессов, то есть Ws=W*s-W-s, 1& lt-5<-^. Если
в стадии участвует третья частица, то скорость 1Л? вычисляется по формулам Г,=Р, Ж. Р,=Т,=1& lt-л<-Ш. где вы, 1 & lt-«<-N5, — эффективности третьих частиц, N1 — число инертных веществ, в5- и сI, NR+1& lt-г<-NR+NI — эффективности и концентрации инертных веществ, соответственно. Значения компонент вектора определяются из схемы химической реакции (10) по соотношениям
V/1 Л=/сЛП, Л"Мс& quot-'-"-. Ъ/5=к-Л г=1ЖЛ где к и к_», 1& lt-«<-^, — константы скоростей прямой и обратной стадий, соответственно. Константы скоростей стадий вычисляются по формулам
к/=А-/ и, сх р {-/•. }/ ЯТ). где Т — температура смеси в реакторе- А/, п! и Е/Я — заданные постоянные. Заметим, что константы скоростей в случае неизотермического реактора постоянными не являются — они зависят от температуры. Однако исторически сначала рассматривался изотермический реактор, и к, 1& lt-г<-^, в настоящее время по-прежнему называют константами.
Стехиометрическая матрица АТ с элементами а, формируется из кинетической схемы (10) по следующему правилу. Номер стадии совпадает с номером столбца, а номер реагента с номером строки матрицы АТ. Если x¦l выступает как исходный реагент, то аг/=-аг/, если x¦l — продукт, то аг/=рг/. Если x¦l является одновременно исходным реагентом и продуктом, то аг/=-аг/+рг/. Обычно в элементарной стадии участвует небольшое количество реагентов, то есть стехиометрическая матрица сильно разрежена.
Численное моделирование пиролиз этана. Пиролиз этана в отсутствии кислорода описывается небольшой последовательностью стадий. Механизм пиролиза этана неоднократно обсуждался в литературе. Здесь принята схема реакции, предложенная и исследованная в [10]:
^63 + ^ СН3+С2Н6^С2Н4+С2Н С2Н5^С2Н4+Н Н + С2Н6 ^H2+C2H5, С2Н5+С2Н5 ^С4Н10'
Константы скоростей стадий имеют вид: k^i. 34−10−5, k2=3. 73−102, ?3=3. 69−103, ?4=3. 66−105 и ?5=1. 62−107. Обозначим концентрации реагентов следующим образом: С2=[СНз], C3=[CH4], C4=[c2H5], C5=[C2H4], C6=[H], c7=[H2] и
c8=[C4H10]. Тогда соответствующая система состоит из 8 обыкновенных дифференциальных уравнений вида C'-=ATV, то есть
Начальная концентрация этана Cj=[C2H6] равна 0. 14, для остальных реагентов концентрации равны нулю.
Расчеты осуществлялись с точностью s=10−2. Эффективность алгоритмов интегрирования оценивалась по числу вычислений правой части if задачи (11) на интервале интегрирования. Численное решение осуществлялось на промежутке [0,0. 26] с начальным шагом А=10'-5. Данная задача удовлетворяет «классическому» определению жесткости. В начале интервала интегрирования наблюдается переходный участок (сотые доли секунды), а затем происходит медленное установление. Сравнение эффективности построенных алгоритмов проводилось с известным методом Мерсона (MERSON) [11]. Для всех методов фактическая точность не хуже задаваемой точности. Алгоритму CESCH42 без контроля устойчивости для нахождения решения потребовалось 22 853 вычислений правой части, для алгоритма с контролем устойчивости CESCH42st имеет место if=20 403, для алгоритма переменного порядка и шага CESCH42vp затраты if=2 588, а для метода Мерсона if=25 796.
Заключение. Из результатов расчетов можно сделать следующие выводы. Во-первых, построенный алгоритм интегрирования второго порядка с контролем точности вычислений и устойчивости численной схемы, а также алгоритм переменного порядка и шага можно применять для решения достаточно жестких задач. Во-вторых, по вычислительным затратам алгоритм переменного порядка и шага CESCH42vp эффективнее метода Мерсона почти в 10 раз. Это является следствием контроля устойчивости численной схемы и расчетов с переменным порядком. Представляется, что при достаточно большой размерности задачи (1) метод CESCH42vp может конкурировать с неявными методами на задачах умеренной жесткости, потому что в нем не обращается матрица Якоби. При решении двенадцати тестовых задач [2] и десяти примеров [12] преимущество алгоритма CESCH42vp выше.
СПИСОК ЛИТЕРАТУРЫ
1. Hairer, E., Norsett, S.P., Wanner, G. Solving ordinary differential equations. Stiff and differential-algebraic problems. Berlin: Springer-Verlag, 1987. 528 p.
2. Hairer, E., Wanner, G. Solving ordinary differential equations. Non stiff problems. Berlin: Springer-Verlag, 1991. 601 p.
3. Novikov, E.A. Explicit methods for stiff systems. Novosibirsk: Nauka, 1997. 197 p.
4. Novikov, E.A., Shornikov, Yu.V. Computer simulation of stiff hybrid systems. Novosibirsk: Novosibirsk State Technikal University Publ. 2013. 451 p.
5. Новиков Е. А., Захаров А. А. Согласование областей устойчивости в явном трехстадийном методе типа Рунге-Кутта // Вестник Тюменского государственного университета. 2011. № 7. Серия «Физико-математические науки. Информатика». С. 187−192.
6. Shampine, L.M. Implementation of Rosenbrock methods. ACM Transaction on Mathematical Software. 1982. № 5. Pp. 93−113.
7. Novikov, V.A., Novikov, E.A. Control of the stability of explicit one — step methods of integrating ordinary differential equations. Soviet Math. Dokl. 1984. Vol. 30, № 1. Pp. 211−215.
8. Ceschino, F., Kuntzman, J. Numerical solution of initial value problems. New Jersey: Prentice-Hall, Englewood Clis, 1966. 318 p.
9. Novikov, E.A. Numerical modeling of a modified oregonator by the (2,1)-method for solving stiff problems. Numerical methods and programming. 2010. Vol. 11, № 2. Pp. 123 130.
10. Kulich, D.M., Taylor, J.E. Mathematical simulation of the oxygen ethane reaction. J. Chem. Kinetic. 1975. Vol. 8. Pp. 89−97.
11. Merson, R.H. An operational methods for integration processes. Australia: Proc. of Symp. on Data Processing. 1957. Pp. 329−331.
12. Enright, W.H., Hull, T.E. Comparing numerical methods for the solutions of systems of ODE'-s. BIT. 1975. № 15. Pp. 10−48.
REFERENCES
1. Hairer, E., Norsett, S.P., Wanner, G. Solving ordinary differential equations. Stiff and differential-algebraic problems. Berlin: Springer-Verlag, 1987. 528 p.
2. Hairer, E., Wanner, G. Solving ordinary differential equations. Non stiff problems. Berlin: Springer-Verlag, 1991. 601 p.
3. Novikov, E.A. Explicit methods for stiff systems. Novosibirsk: Nauka, 1997. 197 p. (in Russian)
4. Novikov E.A., Shornikov Yu.V. Computer simulation of stiff hybrid systems. Novosibirsk: NSTU Publishers, 2013. 451 p. (in Russian)
5. Novikov, E.A., Zaharov A.A. Alignment of Stability in the Explicit Three-Stage Runge-Kutta Method. Vestnik Tjumenskogo gosudarstvennogo universiteta — Tyumen State University Herald. 2011. № 7. Pp. 187−192. (in Russian)
6. Shampine, L.M. Implementation of Rosenbrock methods. ACM Transaction on Mathematical Software. 1982. № 5. Pp. 93−113.
7. Novikov, V.A., Novikov, E.A. Control of the stability of explicit one — step methods of integrating ordinary differential equations. Soviet Math. Dokl. 1984. Vol. 30, № 1. Pp. 211−215.
8. Ceschino, F., Kuntzman, J. Numerical solution of initial value problems. New Jersey: Prentice-Hall, Englewood Clis, 1966. 318 p.
9. Novikov, E.A. Numerical modeling of a modified oregonator by the (2,1)-method for solving stiff problems. Numerical methods and programming. 2010. Vol. 11, № 2. Pp. 123 130.
10. Kulich, D.M., Taylor, J.E. Mathematical simulation of the oxygen ethane reaction. J. Chem. Kinetic. 1975. Vol. 8. Pp. 89−97.
11. Merson, R.H. An operational methods for integration processes. Australia: Proc. of Symp. on Data Processing. 1957. Pp. 329−331.
12. Enright, W.H., Hull, T.E. Comparing numerical methods for the solutions of systems of ODE'-s. BIT. 1975. № 15. Pp. 10−48.

ПоказатьСвернуть
Заполнить форму текущей работой