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

МАТЕМАТИКА И ИНФОРМАТ1

ч_________________________

УДК 519.622 Е.А. Новиков, Л.В. Кнауб

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ОРЕГОНАТОРА ДВУХСТАДИЙНЫМ МЕТОДОМ ТИПА РУНГЕ-КУТТА*

Получено неравенство для контроля устойчивости явного двухстадийного метода типа Рунге-Кутта. Разработан алгоритм переменного порядка и шага, в котором наиболее эффективная численная схема выбирается из критерия устойчивости. Приведены результаты моделирования модифицированного орегонатора.

Ключевые слова: методы Рунге-Кутта, оценка глобальной ошибки, контроль точности вычислений, орегонатор, предельный цикл.

Ye.A. Novikov, L.V. Knaub

OREGONATOR NUMERICAL MODELING BY TWO-STAGE RUNGE-KUTTA TECHNIQUE

The inequality for stability control of explicit two-stage Runge-Kutta technique is obtained in the article. The algorithm of variable order and step in which the most efficient numerical scheme is chosen from the stability criterion is developed. Modeling results of modified oregonator are given.

Key words: Runge-Kutta methods, global error appraisal, calculations exactness control, oregonator, limiting

cycle.

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

Здесь построено неравенство для контроля устойчивости явного двухстадийного метода типа Рунге-Кутта. Разработан алгоритм переменного порядка и шага, в котором наиболее эффективная численная схема выбирается из критерия устойчивости. Алгоритм интегрирования применяется для численного моделирования орегонатора, дающего сложный предельный цикл.

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

У = / у У К =у_ г <г<гк (1)

рассмотрим явный двухстадийный метод типа Рунге-Кутта вида

* Работа поддержана грантами РФФИ №08-01-00621 и президента НШ-3431.2008.9.

Уп+ = У„ + Рк + Р_к_ к =И/ уп к = И/ уп + [3

(2)

где у и / - вещественные М-мерные вектор-функции; t - независимая переменная; h - шаг интегрирования; к1 и к2 - стадии метода;

Р1, р2 и в - числовые коэффициенты, определяющие свойства точности и устойчивости (2). В случае неавтономной системы у'/, у), стадии схемы (2) записываются в виде к1=Л/([п,уп), к1=Ь!(Ъ+№, уп+в^). Ниже для сокращения выкладок будем рассматривать задачу (1). Однако построенные методы можно применять для решения неавтономных задач.

Получим соотношения на коэффициенты метода (2) второго порядка точности. Для этого разложим стадии к\ и кг в ряды Тейлора по степеням /7 до членов с № включительно и подставим в первую формулу (2). При

условии Уп= У (п , сравнивая ряды для точного у([л+1) и приближенного ул+1 решений до членов с Ь2,

запишем условия второго порядка точности схемы (2), которые имеют вид р1+р2=1 и вР2=0,5. При данных соотношениях локальную ошибку 5п схемы (2) можно записать в виде

Построим неравенство для контроля точности вычислений метода второго порядка. Для этого рассмотрим вспомогательную схему уп+1= уп+к1 первого порядка точности. С помощью идеи вложенных методов оценку аналога глобальной ошибки £п,2 метода второго порядка можно вычислить по формуле £п,2=уп+1,1-уп+1, 1 =Р2(к2-). Для повышения надежности данной оценки выберем 0=1. Тогда стадия к вычисляется в точке ^, а к2 в точке ^+1. Как показывают расчеты, использование информации в крайних точках шага приводит к более надежному неравенству для контроля точности вычислений. При 0=1 коэффициенты метода второго порядка определяются однозначно Р1=Р2=0,5, а локальная ошибка 5п и неравенство для контроля точности имеют соответственно вид

где ||* - некоторая норма в

£ - требуемая точность интегрирования.

Теперь построим неравенство для контроля устойчивости (2) предложенным в [2] способом. Для этого рассмотрим вспомогательную стадию кз=Ь/(уп+1). Заметим, что кз совпадает со стадией к, которая применяется на следующем шаге интегрирования и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем стадии к1, к2 и кз применительно к задаче у'=Ау, где А есть матрица с постоянными коэффициентами. В результате получим к1=Хуп, к2=(Х+Х2)уп, к1=(Х+Х2+0,5Х3)уп, где Х=ЬА. Легко видеть, что к2-к1=Хул и 2(кз-к2)=Х3уп. Тогда согласно [3], оценку максимального собственного числа Vn,2=hAmax матрицы Якоби системы (1) можно вычислить по формуле Vn,2=2max1<i<N(|kiз-ki2|l| к2-к'1|). Интервал устойчивости схемы (2) второго порядка точности приблизительно равен двум. Поэтому для ее контроля устойчивости можно применять неравенство Уп,2<2.

Полученная оценка является грубой, потому что вовсе не обязательно максимальное собственное число сильно отделено от остальных, в степенном методе применяется мало итераций, дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг Ьп+1 будем вычислять следующим образом. Новый шаг Ьас по точности определим по формуле Ьас=фЬп, где Ьп есть последний успешный шаг интегрирования, а ф, учитывая соотношение к2-к1=0(Ь2п), задается уравнением д21||к2-к1||=£. Шаг № по устойчивости зададим формулой ^^Ьп, где ф, учитывая соотношение Уп,2=0(Ьп), определяется из равенства ^,2=2. Тогда прогнозируемый шаг вычисляется по формуле Ьп+1=тах[Ьп, тт(Ьзс, Ь^)].

Для численного решения задачи (1) рассмотрим схему вида:

* =-[

+

5п= -(1/12)Л3/ "/ 2+(1/6)Л¥ '2/+0(Л4), 0,5|| к2-к1||<£,

у,» =Уп+ГА + ГА к. - ¥ уп к. =¥ у„+ К

(3)

Заметим, что при Г2=Г1=0,5 численная формула (3) имеет второй порядок точности. Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим (3) для решения скалярного тестового уравнения у'=Ку% Ре(Л)<0. Получим уп+1=0(х)уп, где функция устойчивости 0(х) имеет вид 0(Х)=1+(Г2+Г1)Х+Г2Х2, Х=Ь\. Требование первого порядка точности приводит к соотношению Г2+Г1 = 1, которое ниже будем считать выполненным. Теперь выберем Г2 таким образом, чтобы метод (3) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Т2(г)=212-1 на промежутке [-1,1]. Проведем замену переменных, полагая г=1-2х/у. Получим Т2(х)=1 -8х//+8х2//2, при этом отрезок [/,0] отображается на [-1,1]. Нетрудно показать, что среди всех многочленов вида Р2(х)=1+х+с2х2 для Т2(х) неравенство |Тг(х)|<1 выполняется на максимальном интервале [/,0] при /=-8. Потребуем совпадения коэффициентов О(х) и Т2(х) при /=-8. Это приводит к соотношениям Г2+Г1=1 и Г2=1/8. В результате имеем коэффициенты П=7/8 и /2=1/8 метода порядка точности с максимальным интервалом устойчивости, локальная ошибка 5„ которого имеет ВИД 5 = + . Для контроля точности численной формулы первого порядка

будем использовать оценку локальной ошибки. Учитывая, что к2 - = + и вид локальной

ошибки, неравенство для контроля точности записывается в виде 3||к2-к1||<8£, где £ есть требуемая точность расчетов.

Построим неравенство для контроля устойчивости метода первого порядка. Для этого рассмотрим вспомогательную стадию кз=Ь/(уп+1). Заметим, что кз совпадает со стадией к1, которая применяется на следующем шаге интегрирования и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем стадии к1, кг и кз применительно к задаче у'=Ау. В результате получим к1=Хуп, к2=(Х+Х2)уп, к1=(Х+Х2+0.125Хз)уп, где Х=ЬА. Легко видеть, что к2-к1=Х*уп и 8(кз-к2)=Ху,. Тогда, согласно [3], оценку максимального собственного числа уп,1^Лтах матрицы Якоби можно вычислить по формуле Уп,1=8тах1<,<л/(|к'з-к'2|/| к^к^). Интервал устойчивости схемы (3) первого порядка точности равен 8. Поэтому для ее контроля устойчивости можно применять неравенство Уп,1<8. Шаг выбирается по аналогии с методом второго порядка точности.

На основе построенных методов первого и второго порядков точности легко сформулировать алгоритм переменного порядка и шага. Расчеты всегда начинаются методом второго порядка как более точным. Переход на схему первого порядка осуществляется при нарушении неравенства Уп,2<2. Обратный переход на метод второго порядка происходит в случае выполнения неравенства Уп,1<2. При расчетах по методу первого порядка наряду с точностью контролируется устойчивость.

Расчеты проводились на 1ВМ РС АМоп(ЩХР 2000 с двойной точностью. В конкретных расчетах левая часть неравенства для контроля точности вычислялась по формуле ||к2-к1|| =тахы<ы(|к2-к'1|/|уп|+г), где / - номер компоненты; г - положительный параметр. Если по 1-й компоненте решения выполняется неравенство |уп|<г, то контролируется абсолютная ошибка ге, в противном случае, - относительная ошибка £. В расчетах параметр г выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой. Значение £ полагалось равным 10-2. Это связано с тем, что в построенном алгоритме применяются схемы первого и второго порядков точности, и поэтому проводить расчеты данным методом с более высокой точностью нецелесообразно.

Сравнение эффективности проводилось с тринадцатистадийным методом Рунге-Кутта-Фельберга с контролем ^^788^ [4] и без контроля (РКР78) [5] устойчивости. В [5] получены два комплекта коэффициентов вложенных методов седьмого и восьмого порядков точности, что позволяет оценить ошибку вычислений. Ниже через '&а, 'то и 1/и обозначены соответственно суммарное число шагов интегрирования, число повторных вычислений решения (возвратов) вследствие нарушения требуемой точности расчетов и число вычислений правой части задачи (1).

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

А + УП Х + Р, кх = ),084,

Х + т Р, к2 - • ,

А + ХП IV, к3 = • ,

С + ]¥П X + Zl к4 = ,3-Ю5,

к_ = , к — ■ ,

к_ = • , к_ = >,4-107,

2XU А + Р, к5= • , к_ = • ,

Z^- + >,462Y, к6=), 65,

где к, -5</<6 - константы скоростей прямых (с положительными индексами) и обратных (с отрицательными индексами) стадий. В данной реакции участвуют 7 частиц, имеющих следующие обозначения: A=BrO~3, C=M(n), P=HOBr, W= BrO2, X=HBrO2, Y=Br, Z=M(n+1). В этих обозначениях M(n) - ион металла катализатора, M(n+1) - окисленная форма этого иона. Обозначим концентрации реагентов следующим образом:

C1=[BrO-a], C2=[Br], сз= [M(n)], C4= [HBrO2],

C5=[HOBr], C6= [BrO2], C7=[M(n+1)].

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

1

C' = AVT + ^ .

в

где C=(c1, C2, ..., cj)T- вектор концентраций реагентов; A - стехиометрическая матрица; V=(v1, V2, ..., V6)T - вектор скоростей стадий; Cp=(cp1, Cp2, ..., Cp7)T- вектор концентраций реагентов на входе в реактор; 0 - время пребывания смеси в реакторе; 0=rlu, r - объем реактора; u - объемная скорость течения смеси через реактор. Соответст-

вующая система дифференциальных уравнений имеет вид [5]:

с = —v - v + v + с —с О,

с = -V - V + V + ср_ - с в,

c'=-v+v + cp-c в,

с = v - v_ - v + v - v + ср - с О, (4)

c'=v+ v_+v + с -с О,

С = V -V + ср -с в,

с'' = V -V + ср -с О,

где скорости V1, V2, .„, V6 стадий определяются по формулам:

v, = ксс—к с с , v2= к с с -к с ,

v3 = к с с -к с , v4=k с с —к с с ,

v5= к с -к с с , v6= к с .

Интегрирование системы (4) проводилось на промежутке (0,1000) с начальным шагом 10-5. Концентрации реагентов на входе в реактор принимают значения Cp1=0,14, Cp2=0,151-10-5, Cp3=0,125-10-3,

Cp4=Cp5=Cp6=Cp7=0, причем 0=125,5. Начальные значения концентраций реагентов следующие:

C1=0,1387, C2=0,1534•10-6, C3=0,117610-3, C4=0,3165•10-7,

C5=0,1956 10-3, C6=0,581410-6, C7=0,63110-5.

Алгоритмом переменного порядка и шага RK2PP решение задачи (4) удалось вычислить с затратами /sa=113350, /wo=1878, /fu=228579. Для метода второго порядка RK2 без контроля устойчивости затраты /sa=352636, /wo=117318, ifu=822591, а для алгоритма RK2ST с контролем устойчивости - /sa=309332,

wo=1791, /fu=620456. Вычислительные затраты для метода Фельберга следующие: для RKF78 - /sa=140134, /Wo=140089, /fu=3502810; для RKF78ST - /sa=138912, /wo=46378, ifu=2362392.

Рассмотрим еще простейшую моделью орегонатора с периодическим решением [6], которая в последнее время активно используется для тестирования методов интегрирования:

У' = У. ~ У У. + У. ~ ' У ■

/= -у.-уу.+у 27 ’

у[= У.-у. > (5)

te ,\ = ,J1(°)= ■ ^2(°)= ■ ^з(°)= ■

Алгоритмом RK2PP решение задачи (5) вычислено с затратами /sa=1047298, /Wo=1993, /fu=2096590. Для метода второго порядка RK2 без контроля устойчивости затраты /sa=3735479, /wo=1210999, fu=8681958, а для алгоритма RK2ST с контролем устойчивости - /sa=3913734, /wo=1890, /fu=7829359. Вычислительные затраты для метода Фельберга следующие: для RKF78 - /sa=1 478475, /wo=1477617, /fu=38429196; для RKF7ST- /sa=1507475, /wo=21158, /fu=19872229.

Из сравнения результатов расчетов следует, что контроль устойчивости всегда приводит к повышению эффективности. Это является следствием устранения возвратов (повторных вычислений решения), возникающих из-за неустойчивости численной формулы. В конце интервала интегрирования фактическая точность вычислений для алгоритмов без контроля устойчивости примерно на порядок лучше задаваемой, а для алгоритмов с контролем примерно на два порядка. Такая же тенденция сохраняется при интегрировании всех рассмотренных задач и, в частности, примеров [6-7].

Выводы

Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f. Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет малости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого и с таким шагом осуществляется интегрирование до тех пор, пока не нарушается неравенство для контроля точности. Как правило, число таких шагов невелико. Однако величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторятся многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый.

Литература

1. Shamp/ne, L.M. Implementation of Rosenbrock methods / L.M. Shamp/ne // ACM Transaction on Mathematical Software. - 1982. - №5. - Vol. 8. - P. 93-113.

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

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

4. Fehlberg, E. Classical fifth-, sixth-, seventh- and eighth order Runge-Kutta formulas with step size control / E. Fehlberg// Computing. - 1969. - Vol. 4. - P. 93-106.

5. Новиков, Е.А. Автоматизированный программный комплекс моделирования кинетики сложных химических реакций / Е.А. Новиков, ЮА Шитов, ВИ. Бабушок - Красноярск: ВЦ СО АН СССР, 1990. - 47c.

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

7. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи / Э. Хайрер, С. Нер-сетт, Г. Ваннер. - М.: Мир, 1990. - 512 с.

УДК 539.3 АД. Матвеев

ПРОЦЕДУРА ПОСТРОЕНИЯ ДВУСТОРОННИХ ОЦЕНОК ПОГРЕШНОСТЕЙ КОНЕЧНОЭЛЕМЕНТНЫХ

РЕШЕНИЙ ПЛОСКОЙ ЗАДАЧИ УПРУГОСТИ

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

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

A.D. Matveev TWO-SIDED VALUES CONSTRUCTION PROCESS OF FINITE ELEMENT SOLUTIONS ERRORS OF PLANE ELASTICITY TASK

The numerical process of two-s/ded values construct/on for the relat/ve errors of plane task solut/on of elast/c/ty theory bu/lt us/ng the f/n/te element techn/que by means of the regular part/t/on/ng /s proposed /n the art/cle. Implementat/on of the process /s reduced to the bu/ld/ng of the upper and lower values funct/ons for the relat/ve errors of elast/c gr/d solut/ons of the spec/f/c class of plane elast/c/ty task.

Key words: upper and lower values funct/ons, relat/ve error, homogeneous bod/es, plane elast/c/ty task, f/n/te element techn/que.

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

5 =и0- ||< :0 = (а)

где и0 - точное решение;

Ъ - шаг регулярной узловой сетки дискретной задачи;

и-целое, п> ;

С - константа, которая не зависит от И.

Оценки типа (а) играют важную роль в теоретических исследованиях численных методов [3-4]. Однако применение этих оценок в расчетах связано с трудностями, которые состоят в следующем.

Во-первых, в большинстве случаев С зависит от точного решения или ее частных производных. Поскольку точное решение неизвестно, то не представляется возможным определить значение С .