CC BY
56
M. I. Antipin, I. N. Gusev
MATHEMATICAL MODEL OF THE ECRANOPLAN PARAMETERS CHOICE AT THE STAGE OF THE TECHNICAL OFFER
The mathematical model of the rational parameters choice is constructed at outline designing by a method ofresearch of space of parameters. Values of aerodynamic factorscy=f(a, h) andmz = f(a, h), relative coordinate of aerodynamic focus x= f(a, h), distribution of aerodynamic loading on a bearing surface for three aerodynamic schemes of bearing surfaces plane, «duck», «hybrid» are received numerically. Functional dependences cy = f(a, h), mz = f(a, h), x^ = f(a, h) are constructed.
Кeywords: mathematical model, ground effect, ekranoplan, center of gravity, reversed delta wing.
УДК 519.622
Е. А. Новиков, Л. В. Кнауб
ПРИМЕНЕНИЕ ЯВНОГО ТРЕХСТАДИЙНОГО МЕТОДА ТИПА РУНГЕ-КУТТА ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ЗАДАЧ ХИМИЧЕСКОЙ КИНЕТИКИ1
Получены коэффициенты явного трехстадийного метода типа Рунге-Кутта. Построены неравенства для контроля точности вычислений и устойчивости численной схемы. Приведены результаты расчетов двух моделей орегонатора, подтверждающие повышение эффективности за счет дополнительного контроля устойчивости.
Ключевые слова: методы Рунге-Кутта, контроль устойчивости, орегонатор.
При решении задачи Коши для обыкновенных дифференциальных уравнений средней жесткости и большой размерности в ряде случаев можно применять явные методы. Они не требуют обращения матрицы Якоби и поэтому могут быть эффективнее неявных численных схем. Однако для эффективного использования явных методов при решении задач средней жесткости требуется контролировать не только точность вычислений, но и устойчивость численной схемы. В противном случае, на участке установления вследствие противоречивости требований точности и устойчивости шаг интегрирования раскачивается. В лучшем случае это приводит к большому количеству повторных вычислений решения, а шаг выбирается значительно меньше максимально допустимого. В настоящее время можно выделить два подхода к контролю устойчивости [1; 2]. Первый связан с оценкой максимального собственного числа матрицы Якоби / через ее норму с последующим контролем неравенства А/Ц < Б [1], где h - шаг интегрирования, а положительная постоянная Б зависит от размера области устойчивости. Ясно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к значительному увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числа Хшах матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства А|1| < Б [2]. Такая оценка фактически не приводит к увеличению
затрат [3]. Здесь построено неравенство для контроля устойчивости явного трехстадийного метода типа Рунге-Кутта. Алгоритм интегрирования применяется для численного моделирования орегонатора, дающего сложный предельный цикл.
Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений
у' = /ОХ у('„) = Уо > 1 < Ч > (1)
рассмотрим явный трехстадийный метод типа Рунге-Кут-та, который имеет вид
Уп+1 = Уп + М + Р2 к2 + РзК К = ¥(Уп); К = ¥ (Уп + Р21к1); (2)
к3 = ¥ (Уп + Рз1К1 + Рз2 К2), где у и / - вещественные п-мерные вектор-функции; / - независимая переменная; А - шаг интегрирования; К1, к2 и к3 - стадии метода,р1,р2,р3; Р21, Р31 и Р32 - числовые коэффициенты, определяющие свойства точности и устойчивости. В случае неавтономной системы у = /X у), у('0) = у0, /0 < / < /к, схема (2) записывается следующим образом:
Уп+1 = Уп + Р1к1 + Р 2 к2 + Р3к3; к1 = ¥ ({п > Уп к2 = ¥ Уп + Р2А Уп +Ь21к1); (3)
к3 = ¥ Уп + [Ь31 + Р32 ]А, Уп + Ь31к1 + Р32 к2 ).
Ниже для сокращения выкладок будем рассматривать задачу (1). Однако построенные методы можно применять для решения неавтономных задач.
1 Работа поддержана грантами РФФИ .№08-01-00621 и Президента НШ-3431.2008.9.
Получим соотношения на коэффициенты метода (2) третьего порядка точности. Для этого разложим стадии к1, &2 и к3 в ряды Тейлора и подставим в первую формулу трехстадийного метода (2):
Уп+1 = Уп + (Л + Р 2 + Рз)¥п +
+№21 р2 + (Рзі +Рзі)рз]^2 / 'п/п +
+Л3[Р2,Р32Рз /'/ + 0,5(Р2іР2 +
+(Рзі +Рз2)2 Рз)/\/1] + (4)
+Л4[0,5Р2іРз2 Рз /' п// +
1 +Р21 (Р31 + Р32 )Р32Р3 Г п/'п/П +
+ 6(Р21 р + (Р31 +Р32 )3 Р3)/'"к/3] + 0 (^5).
6
Сравнивая ряды для приближенного уп + 1 и точного у(/и + 1) решений до членов с включительно при усло-
вии уп=у(/п), запишем условия третьего порядка точности схемы (2), которые имеют вид
Р1 + Р2 + Рз = 1; Р21Р2 + (Р31 +Р32 )Рз = 0,5
Р21Р2 + (Р31 + Р32 )2 Р3 = 1/3;
Р21Р32 Рз = 1/6. (5)
Локальную ошибку 5п схемы (2) можно вычислить по формуле 5п=у(ґп + 1) - уп + 1. Учитывая представления у(ґп + 1) и уп + 1 в виде рядов Тейлора, получим
§п = ^ /'Ъ/ + [^-1Р21Р32Рз]///2 +
+[О Р21 (Р31 +Р32 )Р32Р3]///2 +
О
+[2г-6Р21Р2 -6(Р31 +Р32)3Рз]/П + 0(¥).
24 6 6
В нелинейной системе (5) два свободных коэффициента. Исследуем два варианта.
Вариант 1. Минимизируем локальную ошибку выражения (6). Для этого, учитывая его вид, вместо уравнения (5) рассмотрим следующую расширенную нелинейную систему
Р1 + Р2 + Рз = 1;
Р21Р2 + (Р31 +Р32 ) Рз = 1/2;
Р21Р2 + (Р31 +Р32 )2 Р3 = 1/3;
Р21Р32 Рз = 1/6; (7)
Р21Р32 Рз = 1/12;
Р 21 (Р31 +Р32 )Рз2 Рз = 1/8.
При 1,5 Р21 = Р31 + Р32 два последних уравнения (7) совпадают. Из четвертого и пятого соотношений (7) имеем Р21 = 0,5. Из второго и третьего равенств получим р2 = 1/3 и рз = 4/9. Из первого уравнения (7) - р1 = 2/9, а из четвертого - Р32 = 3/4. Наконец, из соотношения Р31 + Р32 = 3/4 получим Р31 = 0. В результате коэффициенты метода (2) с минимальной локальной ошибкой можно записать в виде
р1 = 2/9; р2 = 1/3; р3 = 4/9;
Р21 = 1/2; Р31 = 0; Р32 = 3/4. (8)
При данных коэффициентах локальную ошибку 5 схемы (2) можно представить следующим образом 5п = (/г4/24)х
х/3/- (к4/288)/ ' '/ = 0(к5). При использовании (2) с наборами коэффициентов (8) ни одна стадия не вычисляется в точке 1п +1. При быстром изменении решения это может приводить к понижению точности расчетов.
Вариант 2. Положим Р21 = 0,5 и Р31 + Р32 = 1. Тогда на каждом шаге к,, к. и к, вычисляются в точках /, / + 0,5к и
1’ 2 3 п п ’
I + к, соответственно. В этом случае условия третьего порядка записываются в виде
р1+р2+Рз = 1;
0,5 р2 + р3 = 0,5;
0,25 р2 + р3 = 1/3;
Р32 р3 = 1/3. (9)
Из второго и третьего уравнений данной системы имеем р2 = 2/3 и р3 = 1/6. Из первого и последнего - р1 = 1/6 и Р32 = 2. Из равенства Р31 + Р32 = 1 следует Р31 = -1. В результате коэффициенты метода (2) имеют следующий вид
р1 = 1/6; р2 = 2/3; р3 = 1/6;
Р 21 = 0,5; Р3, =-1; Р32 = 2. (10)
При данных соотношениях локальную ошибку 5п схемы (2) можно записать следующим образом 5 = (к4/24)х
х /3/-/ '// - (1/3)/ ' '/] = 0(к5).
Построим неравенство для контроля точности вычислений метода третьего порядка. Для этого рассмотрим вспомогательную схему
Уп+1,1 = Уп + Г1к1 + Г2к2 > (11)
(6) где к1 и к2 определены в уравнениях (2). Необходимо, чтобы данный метод имел второй порядок точности. Сравнивая ряды для у(/п + 1) и уп + 11 видим, что требование второго порядка будет выполнено, если г + г2 = 1, Рг2=0,5. Отсюда получим г2 = 0,5/Р, г1 = 1-г2, где значение в21 определено в (8) или (10). С помощью идеи вложенных методов оценку ошибки еп3 метода третьего порядка можно одштъ по формуле 8п,3 = уп + ! - уп + 1Д = (р 1 - г1)к1 + р - г2) к2 + р3к3. Тогда неравенство для контроля точности вычислений имеет вид ||(р1-г1)к1 + (р2-г2)к2 + р3к3|| < е, где ||-|| - некоторая норма в Яп, е - требуемая точность интегрирования. В конкретных расчетах применялся метод (2) с коэффициентами (10), как более надежный. Исходя из этого неравенство для контроля точности вычислим следующим образом: ||к1 - 2к2 + к3|| < 6е.
Теперь построим неравенство для контроля устойчивости (2) предложенным в работах [2; 3] способом. Запишем стадии к1, к2 и к3 применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами. В результате получим
кх = Хуп; к2 = (X + Р21X2)уп;
к3 = [ X + (Р3, +Р32) X2 + Р21Р32 X3] уп, гдеX = кА. Найдем коэффициенты йр ^2 и ^3 из условия выполнения равенства
йхкх + ^2^2 + ^3к3 = X Уп.
Данное требование будет выполнено, если
й1 = (Р31 + Р32 Р21 )/ й;
й2 = -(Р31 +Р32)/ й ; й3 =Р21 / й, где й = Р221Р32. Также видно, что
Р-!(к2 - К) = X2y„ ■
Тогда согласно работе Е. А. Новикова [3] оценку максимального собственного числа v = hl матрицы Яко-
n,3 max А
би системы (!) можно вычислить по формуле
vn 3 = b2! max(d!k'! + d2к'2 + d3k'3 |/| К2 - k[ |). (!2)
’ !<'<N
Интервал устойчивости численной схемы (2) приблизительно равен 2,5. Поэтому для ее контроля устойчивости можно применять неравенство vn 3 < 2,5.
Полученная оценка выражения (!2) является грубой, потому что вовсе не обязательно максимальное собственное число сильно отделено от остальных, в степенном методе применяется мало итераций, дополнительные искажения вносит нелинейность задачи (!). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг hn + ! будет вычисляться следующим образом: новый шаг hac по точности определим по формуле hac = q!hn, где hn - последний успешный шаг интегрирования; q учитывая соотношение en3 = O(h3n), задается уравнением q3!|| 3|| = e. Шаг hst по устойчивости запишем формулой hst = q2hn, где q2, учитывая равенство vn3 = O(hn), определяется из уравнения q2vn3 = 2,5. В результате, шаг hn + ! вычисляется по следующей формуле:
hn+i = max[hn, min(hac, hst)].
Данное выражение позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Собственно говоря, именно наличие данного участка существенно ограничивает возможности применения явных методов для решения жестких задач.
Расчеты проводились на Intel(R) Core 2 Quad CPU с двойной точностью. В конкретных расчетах левая часть неравенства для контроля точности вычислялась по формуле ||k2 - к!|| = max!<'<N(|k2 - kj/yy + r), где ' - номер компоненты; r - положительный параметр. Если по '-й компоненте решения выполняется неравенство \у | < r, то контролируется абсолютная ошибка ге, в противном случае - относительная ошибка e. В расчетах параметр r выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой. Значение e полагалось равным !0-2.
Сравнение эффективности проводилось с тринадцатистадийным методом Рунге-Кутта-Фельберга с контролем (RKF78ST) [4] и без контроля (RKF78) [5] устойчивости. В работе E. Fehlberg [5] получены два комплекта коэффициентов вложенных методов седьмого и восьмого порядков точности, что позволяет оценить ошибку вычислений. Ниже через 'sa, 'wo и 'fu обозначены, соответственно, суммарное число шагов интегрирования, число повторных вычислений решения (возвратов) вследствие нарушения требуемой точности расчетов и число вычислений правой части задачи (!).
Сравнение эффективности методов проводилось на модели модифицированного орегонатора, дающего сложный предельный цикл. Соответствующая реакция состоит из шести стадий следующего вида:
A + Y € X + P ; К = 0,084; к-! = !04;
: 2W; k3 = 2 -103; k-3 = 2-107;
: X + Z; k4 = 1,3-105; k-4 = 2,4-107 >A + P; k5 = 4 -107; k-5 = 4-10-11;
X + Y € 2P ; k2 = 4-10°
; k-2 = 5
10-
А + X :
с+w<
2 Х<
г ® С + 0,4627; к6 = 0,65, где к. -5 < і < 6 - константы скоростей прямых (с положительными индексами) и обратных (с отрицательными индексами) стадий. В данной реакции участвуют 7 частиц, имеющие следующие обозначения: А = БгО3, С=М(п), Р=НОБг, W=БгО2, Х= НБгО2, ¥=Бг, г=М(п + 1). В этих обозначениях М(п) - ион металла катализатора, М(п + 1) - окисленная форма этого иона. Обозначим концентрации реагентов следующим образом:
С = [БгО-3], с = [Бг], Сз = [М(п)], С4 = [НБгО^,
с. = [НОБг], с6 = [БгО2], с7 = [М(п + 1)].
Данная реакция протекает в изотермическом реакторе постоянного объема с обменом вещества, то есть соответствующая система дифференциальных уравнений состоит из семи уравнений вида (например, [6; 7])
С = АУТ + 1(СР - С),
0 р
где С = (с1, с2, ..., с7)Т- вектор концентраций реагентов; А - стехиометрическая матрица; V = (ур у2, ..., у6)т - вектор скоростей стадий; Ср = (с , ср2, ..., ср7)Т- вектор концентраций реагентов на входе в реактор; © - время пребывания смеси в реакторе, ©= г/и; г - объем реактора; и - объемная скорость течения смеси через реактор.
Соответствующая система дифференциальных уравнений имеет вид
с1 = -^ - V + V + (ср1- с1)/0;
с2 =-у1- V + °,462п + (ср 2- с2)/0;
с3 = V + V + (срз-сз)/0;
с4 = V - V - V + ^ - 2^5 + (ср4 - с4)/0 ; (13)
с5 = V! + 2V2 + V + (ср5 -с5)/0 ;
с6 = 2^- V + (ср 6- с6)/0;
с7 = V - V + (ср7 - с7)/0 , где скорости У1, У2, ., у6 стадий определяются по следующим формулам:
V = к1с1с2 - к-1с4 с5; у2 = к2 с2 с4 - к-2 с52;
У3 = кзс1с4 - к-3с6 ; У4 = к4сзс6 - к-4с4с7 ;
V. = к5с42 - к-5с1с5; У6 = к6с7.
Интегрирование системы (13) проводилось на промежутке [0,1000] с начальным шагом 10-5. Концентрации реагентов на входе в реактор принимают значения с , = 0,14, с 2=0,151 • 10-5, с 3=0,125 • 10~3, с 4 = с . = с,=с 7=0,
р1 7 7 р2 7 7 р3 7 7 р4 р5 р6 р! 7
причем © = 125,5. Начальные значения концентраций реагентов следующие:
с1 = 0,138 7; с2 = 0,153 4 • 10-6; с3 = 0,117 6 • 103; с4 = 0,316 5 • 107;
3 ’ ’ 4 ’ ’
с. = 0,195 6 • 103; с6 = 0,581 4 • 10-6; с7 = 0,631 • 105.
5 7 7 6 7 7^7
При решении задачи (13) методом третьего порядка ИК3 без контроля устойчивости затраты іі'а = 274 247;
В4
1м>о = 71 808; /и = 966 357, а для алгоритма КК38Т с контролем устойчивости - Ш = 233 770; то = 3 517; 1/и = 708 344. Вычислительные затраты для метода Фельберга следующие: для ЯКБ78 - 1яа = 140 134; 1то = 140 089; 1/и = 3 502 810; для ККР788Т - 1за = 138 912; 1то = 46 378; 1/и = 2 362 392.
Рассмотрим простейшую моделью орегонатора с периодическим решением [8], которая в последнее время активно используется для тестирования методов интегрирования
у' = 77,27(у2 - у,у2 + у, - 8,375 • 10-6 у,2);
у2 =(-у 2- у у 2 + у3)/77,27;
у3 = 0,161(у, - у3); (14)
г е [0,300]; к0 = 10-3; у1(0) = 4 ;
у2(0) = 1,1; у3(0) = 4.
При решении задачи (14) методом третьего порядка КК3 без контроля устойчивости затраты составляют isa = 2 903 698; то = 769 326; 1/и = 10 249 746, а для алгоритма КК38Т с контролем устойчивости - 1$а = 2 966 743; то = 7 764; 1/и = 8 915 757. Вычислительные затраты для метода Фельберга следующие: для ЯКР78 - isa = 1 478 475; 1то = 1 477 617; 1/и = 38 429 196; для ЯКБ78Т -isa = 1 507 475; 1то = 21 158; 1/и = 19 872 229.
Из сравнения результатов расчетов следует, что контроль устойчивости всегда приводит к повышению эффективности. Это является следствием устранения возвратов (повторных вычислений решения), возникающих из-за неустойчивости численной формулы. В конце интервала интегрирования фактическая точность вычислений для алгоритмов без контроля устойчивости примерно на порядок лучше задаваемой, а для алгоритмов с контролем -примерно на два порядка. Такая же тенденция сохраняется при интегрировании всех рассмотренных задач и, в частности, примеров [8; 9].
Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции /. Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет ма-
лости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого и с таким шагом осуществляется интегрирование до тех пор, пока не нарушается неравенство для контроля точности. Как правило, число таких шагов невелико. Однако величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторяться многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый.
Библиографический список
!. Shampine, L. M. Implementation ofRosenbrock method / L. M. Shampine // ACM Transaction on Mathematical Software. !982. Vol. 8. № 5. P. 93-!!3.
2. Новиков, Е. А. Контроль устойчивости явных методов интегрирования обыкновенных дифференциальных уравнений / Е. А. Новиков, В. А. Новиков // ДАН СССР. !984. Т. 277. № 5. С. !058-!062.
3. Новиков, Е. А. Явные методы для жестких систем / Е. А. Новиков. Новосибирск : Наука, !997.
4. Новиков, Е. А. Контроль устойчивости метода Фель-берга седьмого порядка точности / Е. А. Новиков, Ю. В. Шорников // Вычислительные технологии. 2006. Т. !!. № 4. С. 65-72.
5. Fehlberg, E. Classical fifth-, sixth-, seventh- and eighth order Runge-Kutta formulas with step size control / E. Fehlberg // Computing. !969. Vol. 4. Р 93-!06.
6. Бабкин, В. С. Автоматизированный банк кинетической информации (общее описание) / В. С. Бабкин, В. И. Бабушок, Ю. П. Дробышев, Ю. Н. Молин, Е. А. Новиков, Г. И. Скубневская // Препринт № 704. Новосибирск : ВЦ СО АН СССР. !987.
7. Бабкин, В. С. Разработка банка кинетической информации / В. С. Бабкин, В. И. Бабушок, Ю. Н. Молин, Е. А. Новиков, Г. И. Скубневская // Гос. служба стандартных справочных данных. М. : Изд-во стандартов, !989. Бюлл. 20. С. 8-!0.
8. Enright, W. H. Comparing numerical methods for the solutions of systems of ODE’s / W. H. Enright, T. E. Hull // BIT. !975. № !5. Р !0-48.
9. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. М. : Мир, !999.
E. A. Novikov, L. V Knaub
APPLICATION OF EXPLICIT THREE-STAGE RUNGE-KUTTA METHOD FOR NUMERICAL MODELING OF CHEMICAL KINETICS PROBLEMS
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. Numerical results of oregonator two models with an additional stability control demonstrate an efficiency increase.
Keywords: Runge-Kutta methods, stability control, oregonator
