Применение явного трехстадийного метода типа Рунге-Кутта для численного моделирования задач химической кинетики

При решении задачи Коши для обыкновенных дифференциальных уравнений средней жесткости и большой размерности в ряде случаев можно применять явные методы. Они не требуют обращения матрицы Якоби и поэтому могут быть эффективнее неявных численных схем. Однако для эффективного использования явных методов при решении задач средней жесткости требуется контролировать не только точность вычислений, но и устойчивость численной схемы. В противном случае, на участке установления вследствие противоречивости требований точности и устойчивости шаг интегрирования раскачивается. В лучшем случае это приводит к большому количеству повторных вычислений решения, а шаг выбирается значительно меньше максимально допустимого. В настоящее время можно выделить два подхода к контролю устойчивости [1; 2]. Первый связан с оценкой максимального собственного числа матрицы Якоби fy через ее норму с последующим контролем неравенства hIf1|| < D [1], где h - шаг интегрирования, а положительная постоянная D зависит от размера области устойчивости. Ясно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к значительному увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числа X матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства h|Xmax| < D [2]. Такая оценка фактически не приводит к увеличению затрат [3]. Здесь построено неравенство для контроля устойчивости явного трехстадийного метода типа Рунге-Кутта. Алгоритм интегрирования применяется для численного моделирования орегонатора, дающего сложный предельный цикл.

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

У ' = f ( У ) , У ( t 0 ) = У 0 , t 0 <  t <  t k ,         (1)

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

У п + 1 = У п ⁺ P 1^k 1 ⁺ P 2 ^k 2 ⁺ P 3 ^k 3 ;

К = hf ( У п ); k 2 = hf ( У п +в 21 k i );           (2)

k3 = hf (Уп +вз1 ki +вз2 k 2), где y и f – вещественные n-мерные вектор-функции; t – независимая переменная; h – шаг интегрирования; k 1, k2 и k3 - стадии метода, p 1, p2,p3; P21, Р31 и P32 - числовые коэффициенты, определяющие свойства точности и устойчивости. В случае неавтономной системы у' = ft, y), У(10) = У0, 10 < t < tk, схема (2) записывается следующим образом:

У п + 1 = У п ⁺ P 1 ^k 1 ⁺ P 2 ^k 2 ⁺ P 3 ^k 3 ;

k 1 = hf ( t n , У п );, k 2 = hf ( t n +P 21 h , У п +P 21 k 1 );      (3)

k 3 = hf ( t n + [ P 31 +P 32 ] h , У п +0 31 k 1 +P 32 k 2 ).

Ниже для сокращения выкладок будем рассматривать задачу (1). Однако построенные методы можно применять для решения неавтономных задач.

Получим соотношения на коэффициенты метода (2) третьего порядка точности. Для этого разложим стадии к 1 , к ₂ и к ₃ в ряды Тейлора и подставим в первую формулу трехстадийного метода (2):

Ут + 1 = Ут ^{+ (} Р\ ⁺ Р 2 ⁺ Р з ) ^hfm ⁺

+ [ р 21 Р 2 + ( Р *1 + Р *1 ) Р 3 ] h ² f' m f m +

+ h ³[ в 21 Р з2 Р 3 f ff n + 0,5( P 21 Р 2 +

+ ( Р 31 +P 32 )² Р з ) f "  m f m ²] +                 (4)

+ h ⁴ [0,5 Р 21 Р з2 Р 3 f' J "  m f m ' +

+ P 21 ( P 31 +P 32 ) P 32 Р 3 f " m f 'm f 2 +

+ 1( P ³ 21 р 2 + ( P 31 +P 32 )³ Р 3 ) f ^m „ Л '] + O ( h 5 ).

Сравнивая ряды для приближенного y_{m +} 1 и точного у ( t_{m +} 1 ) решений до членов с h ³ включительно при условии у = у ( t ), запишем условия третьего порядка точности схемы (2), которые имеют вид

Р 1 + Р 2 + Р 3 = 1; Р 21 Р 2 + ( Р 31 +Р 32 ) Р 3 = 0,5;

Р 21 p 2 + ( Р 31 +Р 32 )² Р 3 = 1/3;

Р 21 Р 32 Р 3 = 1/6.                     (5)

Локальную ошибку 5 _т схемы (2) можно вычислить по формуле 5 _К = у ( t_{m +} ) - y_{m +} 1 . Учитывая представления у ( t_{m +} ) и у_{т +} 1 в виде рядов Тейлора, получим

5 т = h '{^^ f ³f + [^ - 2 Р 21 Р 32 Р 3 ] fff ² +

+ [2 -p 21 ( p 31 +p 32 ) p 32 Р 3 ] f ff ² +            (6)

■\, - ¹ в ³ 21 Р 2 - 7^ 31 + ₽ 32 )³ Р 3 ] f f ³ } + O ( h ⁵ ).

24 6       6

В нелинейной системе (5) два свободных коэффициента. Исследуем два варианта.

Вариант 1. Минимизируем локальную ошибку выражения (6). Для этого, учитывая его вид, вместо уравнения (5) рассмотрим следующую расширенную нелинейную систему

Р 1 + р 2 + р з = 1;

Р 21 Р 2 + ( Р 31 +Р 32 ) Р 3 = 1/2;

Р 21 Р 2 + ( Р 31 +Р 32 )² р 3 = 1/3;

Р 21 Р 32 Р 3 = 1/6;                      (7)

Р 21 Р з2 Р 3 = 1/12;

p 21 ( p 31 +Р з2 ) Р з2 Р 3 = 1/8.

При 1,5 Р21 = Р31+ Р32 два последних уравнения (7) совпадают. Из четвертого и пятого соотношений (7) имеем Р21 = 0,5. Из второго и третьего равенств получимp2= 1/3 иp* = 4/9. Из первого уравнения (7) -p 1 = 2/9, а из четвертого - Р32 = 3/4. Наконец, из соотношения Р31 + Р32 = 3/4 получим Р31 = 0. В результате коэффициенты метода (2) с минимальной локальной ошибкой можно записать в виде р1 = 2/9; p2 = 1/3; p3 = 4/9;

Р 21 = 1/2; Р 31 = 0; Р 32 = 3/4.               (8)

При данных коэффициентах локальную ошибку 5m схемы (2) можно представить следующим образом 5т = (h4/24)х хf3f— (h4/288)f "f = O(h5). При использовании (2) с наборами коэффициентов (8) ни одна стадия не вычисляется в точке tm + 1. При быстром изменении решения это может приводить к понижению точности расчетов.

Вариант 2. Положим Р ₂₁ = 0,5 и Р ₃₁ + Р ₃₂= 1. Тогда на каждом шаге к 1 , к ₂ и к * вычисляются в точках t_m , t_m + 0,5 h и t_m + h , соответственно. В этом случае условия третьего порядка записываются в виде

Р 1 + Р 2 ⁺ Р з = 1;

0,5 p 2 + p * = 0,5;

0,25 p ₂ + p ₃ = 1/ 3 ;

Р 32 Р 3 = 1/3.                      (9)

Из второго и третьего уравнений данной системы име-ем p ₂= 2/3 и p ₃= 1/6. Из первого и последнего - p 1 = 1/6 и Р ₃₂=2. Из равенства Р ₃₁ + Р ₃₂= 1 следует Р ₃₁=-1. В результате коэффициенты метода (2) имеют следующий вид

Р 1 = 1/ 6; p ₂ = 2 / 3 ; p ₃ = 1/ 6 ;

Р 21 = 0,5; Р 31 =- 1; Р 32 = 2.           (10)

При данных соотношениях локальную ошибку 5 _m схемы (2) можно записать следующим образом 5 = ( h ⁴/24) х х [ f³f - f Tf - (1/3) f "f ³] = O( h ⁵).

Построим неравенство для контроля точности вычислений метода третьего порядка. Для этого рассмотрим вспомогательную схему уп+1,1 = уп + r1к1 + r2 к 2,              (11)

где к ₁ и к ₂ определены в уравнениях (2). Необходимо, чтобы данный метод имел второй порядок точности. Сравнивая ряды для у ( t_{m +} 1 ) и у_{п +} 11 видим, что требование второго порядка будет выполнено, если r 1 + r ₂= 1, Р r ₂= 0,5. Отсюда получим r ₂= 0,5/ р , r 1 = 1- r ₂, где значение в₂₁ определено в (8) или (10). С помощью идеи вложенных методов оценку ошибки s _{m 3} метода третьего порядка можно оценить по формуле ^s m ,3 = у_п +1^- у_к +1,1 = (р 1^- r 1 ) ^к 1⁺ ( р 2^- r 2 ) к ₂ + p ₃ к ₃. Тогда неравенство для контроля точности вычислений имеет вид ||( p 1 - r 1 ) к 1 + ( p ₂- r ₂) к ₂+ p ₃ к ₃|| < s , где || - || - некоторая норма в R^m , s - требуемая точность интегрирования. В конкретных расчетах применялся метод (2) с коэффициентами (10), как более надежный. Исходя из этого неравенство для контроля точности вычислим следующим образом: || к 1 - 2 к ₂+ к ₃|| <  6 s .

Теперь построим неравенство для контроля устойчивости (2) предложенным в работах [2; 3] способом. Запишем стадии к 1, к2 и к3 применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами. В результате получим к1 = ХуК; к2 = (X + Р21X2)уп;

к* = [ X + (Р31 +Р32) X2 +Р21Р32 X3] ук , где X = hA. Найдем коэффициенты d 1, d2 и d3 из условия выполнения равенства d 1 к1 + d 2 к 2 + d 3 к3 = X3 уК.

Данное требование будет выполнено, если d 1 = (Р31 +Рз2-Р21)/d;

d 2 = - ( Р 31 + Р 32 ) / d ; d * = P 21 / d , где d = P ² ₂₁ P ₃₂. Также видно, что

β 2 ^- 1 ¹( k 2 - k 1 ) = X ² y n .

Тогда согласно работе Е. А. Новикова [3] оценку максимального собственного числа νn,3= hλmax матрицы Якоби системы (1) можно вычислить по формуле vn3=β21max(d1k1i+d2k2i+d3k3i|/|k2i-k1i|). (12) ,           1≤i≤N

Интервал устойчивости численной схемы (2) приблизительно равен 2,5. Поэтому для ее контроля устойчивости можно применять неравенство ν _{n ,3} ≤ 2,5.

Полученная оценка выражения (12) является грубой, потому что вовсе не обязательно максимальное собственное число сильно отделено от остальных, в степенном методе применяется мало итераций, дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг h_{n + 1} будет вычисляться следующим образом: новый шаг h^ac по точности определим по формуле h^ac = q ₁ h_n , где h_n – последний успешный шаг интегрирования; q ₁, учитывая соотношение ε _{n ,3} = O ( h ³ _n ), задается уравнением q ³₁|| ε _{n ,3}|| = ε . Шаг h^st по устойчивости запишем формулой h^st = q ₂ h_n , где q ₂, учитывая равенство v_{n ,3} = O ( h_n ), определяется из уравнения q ₂ v_{n ,3} = 2,5. В результате, шаг h_{n + 1} вычисляется по следующей формуле: h_{n + 1} = max[ h_n , min( h^ac , h^st )] .

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

Расчеты проводились на Intel(R) Core 2 Quad CPU с двойной точностью. В конкретных расчетах левая часть неравенства для контроля точности вычислялась по формуле ||k ₂ – k ₁ || = max_{1 ≤ i ≤ N} (| kⁱ ₂ – kⁱ ₁|/| yⁱ_n | + r ), где i – номер компоненты; r – положительный параметр. Если по i -й компоненте решения выполняется неравенство | yⁱ_n | <  r , то контролируется абсолютная ошибка r е, в противном случае – относительная ошибка ε . В расчетах параметр r выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой. Значение ε полагалось равным 10^–2 .

Сравнение эффективности проводилось с тринадцатистадийным методом Рунге-Кутта–Фельберга с контролем (RKF78ST) [4] и без контроля (RKF78) [5] устойчивости. В работе E. Fehlberg [5] получены два комплекта коэффициентов вложенных методов седьмого и восьмого порядков точности, что позволяет оценить ошибку вычислений. Ниже через isa , iwo и ifu обозначены, соответственно, суммарное число шагов интегрирования, число повторных вычислений решения (возвратов) вследствие нарушения требуемой точности расчетов и число вычислений правой части задачи (1).

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

A + Y ^ X + P ; k 1 = 0,084; k - 1 = 10⁴;

X + Y ^ 2 P ; k 2 = 4 - 10^s; k _{- 2} = 5 - 10 - 5 ;

A + X ^ 2W ; k = 2 - 10 3 ; k, = 2 - 10 7 ; 3                  - 3

C + W ^ X + Z ; k ₄ = 1,3 - 10 5 ; k -4 = 2,4 - 10 7 ;

2 X ^ A + P ; k 5 = 4 - 10 7 ; k -5 = 4 - 10 ^{- 11};

Z→C+0,462Y;k6=0,65, где ki –5 ≤ i ≤ 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) – окисленная форма этого иона. Обозначим концентрации реагентов следующим образом:

с1= [BrO–3], с2= [Br–], с3= [M(n)], с4= [HBrO2], с5= [HOBr], с6= [BrO2], с7= [M(n + 1)].

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

C′=AVT+1(Cp-C), θ где C = (c1, c2, …, c7)T – вектор концентраций реагентов; A – стехиометрическая матрица; V = (v1, v2, …, v6)T – вектор скоростей стадий; Cp = (cp1, cp2, …, cp7)T – вектор концентраций реагентов на входе в реактор; Θ – время пребывания смеси в реакторе, Θ= r/u; r – объем реактора; u – объемная скорость течения смеси через реактор.

Соответствующая система дифференциальных уравнений имеет вид c1′=-v1-v3+v5+(cp1-c1)/θ;

c ₂ ′=- v ₁ - v ₂ + 0,462 v ₆ + ( c_{p 2} - c ₂)/ θ ;

c ₃ ′=- v ₄ + v ₆ + ( c_{p 3} - c ₃)/ θ ;

c ₄ ′= v ₁ - v ₂ - v ₃ + v ₄ - 2 v ₅ + ( c_{p 4} - c ₄)/ θ ;     (13)

c ₅ ′= v ₁ + 2 v ₂ + v ₅ + ( c_{p 5} - c ₅)/ θ ;

c ₆ ′= 2 v ₃ - v ₄ + ( c_{p 6} - c ₆)/ θ ;

c7′=v4-v6+(cp7-c7)/θ, где скорости v1, v2, …, v6 стадий определяются по следующим формулам:

v ₁ = k ₁ c ₁ c ₂ - k _{- 1} c ₄ c ₅; v ₂ = k ₂ c ₂ c ₄ - k _{- 2} c ₅;

v ₃ = k ₃ c ₁ c ₄ - k _{- 3} c ₆²; v ₄ = k ₄ c ₃ c ₆ - k _{- 4} c ₄ c ₇;

v ₅ = k ₅ c ₄ ² - k _{- 5} c ₁ c ₅; v 6 = k 6 c 7 .

Интегрирование системы (13) проводилось на промежутке [0,1000] с начальным шагом 10^–5. Концентрации реагентов на входе в реактор принимают значения c_{p 1}= 0,14, c_{p 2}= 0,151 ⋅ 10^–5, c_{p 3} = 0,125 ⋅ 10^–3, c_{p 4}= c_{p 5}= c_{p 6}= c_{p 7}= 0, причем Θ = 125,5. Начальные значения концентраций реагентов следующие:

c ₁= 0,138 7; c ₂= 0,153 4 ⋅ 10^–6;

c ₃= 0,117 6 ⋅ 10^–3; c ₄= 0,316 5 ⋅ 10^–7;

c ₅= 0,195 6 ⋅ 10^–3; c ₆= 0,581 4 ⋅ 10^–6; c ₇= 0,631 ⋅ 10^–5.

При решении задачи (13) методом третьего порядка RK3 без контроля устойчивости затраты isa = 274 247;

iwo = 71 808; ifu = 966 357, а для алгоритма RK3ST с контролем устойчивости – isa = 233 770; iwo = 3 517; ifu = 708 344. Вычислительные затраты для метода Фельберга следующие: для RKF78 – isa = 140 134; iwo = 140 089; ifu = 3 502 810; для RKF78ST – isa = 138 912; iwo = 46 378; ifu = 2 362 392.

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

• У ’ = 77,27( y 2 - y y 2 + y - 8,375 ■ 10 ^{- 6}y²);

• У ‘ = ^{( -} У 2 ^- У 1 У 2 ⁺ У з^)/77,2 7;

У ‘= 0,161( У 1 - у з ); (14)

t e [0,300]; h о = 10 ^{- 3}; У 1 (0) = 4;

y 2 (0) = 1,1; y з (0) = 4.

При решении задачи (14) методом третьего порядка RK3 без контроля устойчивости затраты составляют isa = 2 903 698; iwo = 769 326; ifu = 10 249 746, а для алгоритма RK3ST с контролем устойчивости – isa = 2 966 743; iwo = 7 764; ifu = 8 915 757. Вычислительные затраты для метода Фельберга следующие: для RKF78 – isa = 1 478 475; iwo = 1 477 617; ifu = 38 429 196; для RKF7ST – isa = 1 507 475; iwo = 21 158; ifu = 19 872 229.

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

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