Sensitivity (ЭВМ для спецгруппы), страница 2

Жесткость дифференциальных уравнений химической кинетики не позволяет применять для решения прямой задачи стандартные методы численного интегрирования, такие, как метод Рунге-Кутта, Адамса и т. п. Только в начале 70‑х годов Гиром [21‑23] был разработан эффективный алгоритм интегрирования жестких систем, который сразу же стал общепринятым методом решения кинетических уравнений [19, 20, 24‑27].

Метод Гира представляет собой метод прогноза и коррекции, построенный по схеме Нордсика [28, 29]. Способность интегрировать жесткие уравнения обусловлена двумя его особенностями:

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

б) процедура коррекции реализована как решение системы нелинейных уравнений методом Ньютона. Этот подход обеспечивает существенно более высокую устойчивость решений по сравнению с традиционной схемой коррекции, применяемой, например, в методах типа Адамса.

В соответствии со схемой Нордсика, информация о локальном поведении решения в окрестности текущей точки хранится в виде массива промасштабированных производных :

, (6)

где q – текущий порядок метода, h – текущий шаг интегрирования. Строго говоря, здесь только первая производная является точной, поскольку вычисляется непосредственно по уравнениям (1). Вместо истинных значений высших производных используются их конечно-разностные оценки для q последних шагов. При этом гарантируется, что если локальная погрешность решения в точке t_n (после того, как сделан шаг длиной h) не превышает заданной величины , то в любой точке решение в пределах точности  можно получить с помощью интерполяции по формуле

, (7)

где

.

Вообще говоря, решение можно экстраполировать по формуле (7) и за пределы отрезка , но в этом случае ошибка может превышать величину .

Метод Гира всегда стартует как метод первого порядка (q), так что в начальный момент времени t₀ массив Z целиком определяется начальными условиями и видом функций f_i системы (1). Опуская детали алгоритма, очередной шаг интегрирования можно схематически представить следующим образом:

1) Прогноз. Предполагается, что массив Z содержит компоненты решений и их производные для последней полученной точки t_n (n). Текущий порядок метода равен q, предполагаемая величина шага – h. Оценка решений и масштабированных производных в точке получается путем умножения матрицы Z справа на треугольную матрицу биномиальных коэффициентов (матрицу Паскаля) размером (q)(q):

, (8)

.

Тильда над в левой части (8) означает, что полученный результат не является окончательным, а рассматривается как первоначальная оценка (прогноз) решений и их производных в точке t_n+1.

2) Коррекция. Элементы содержат прогнозируемые значения первых производных в точке t_n+1. Значения этих производных можно также получить непосредственно из уравнений (1), подставляя в правые части оценки решений , причем результат, как правило, будет отличаться от , и это различие связано с погрешностью прогноза. Исправленное решение пытаются получить из уравнений

, (9)

где – некоторый коэффициент, зависящий от порядка метода q, а – искомый вектор решения в точке t_n+1. Так как неизвестные величины в правой части уравнений (9) входят (вообще говоря) в нелинейные функции f_i, эту систему решают методом Ньютона, используя прогноз в качестве начального приближения. При этом на каждой итерации поправки к очередному приближению находятся из системы линейных алгебраических уравнений

(10)

где e – вектор искомых поправок,
 – вектор невязок системы (9) при подстановке в нее текущего приближения,
J – матрица Якоби системы (1) с элементами ,
E – единичная матрица.

Если итерации метода Ньютона не сходятся, это означает, что погрешность прогноза слишком велика. В этом случае уменьшают величину шага h и повторяют попытку, начиная с этапа 1.

3) Оценка погрешности решения. После успешной коррекции вычисляют разницу между прогнозом и окончательным решением и на основании этой величины оценивают фактическую погрешность полученного решения. Пользуясь известными значениями производных, одновременно определяют, какой была бы погрешность решения, если бы шаг той же длины был сделан при порядке q или q. Если оцененная погрешность превышает заданный предел , шаг считается неудачным, величина h уменьшается, и вся процедура повторяется, начиная с прогноза. Если же погрешность меньше, чем , шаг считается успешным, и рассматривается возможность увеличения h на следующем шаге. Это делается путем сравнения трех оценок погрешностей для порядков q, q, q и предельно допустимой погрешности . Будет выбран наиболее выгодный порядок метода, позволяющий максимально увеличить шаг, а, следовательно, ускорить интегрирование при сохранении заданной точности решений.

Отметим некоторые особенности метода Гира, которые необходимо учитывать при его практическом применении:

1) Как уже было сказано, процесс интегрирования начинается при порядке q. Это обычно приводит к весьма малым значениям шага. В дальнейшем порядок меняется, и шаг постепенно растет, увеличиваясь на несколько порядков.

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

I.3. Методы расчета коэффициентов чувствительности.

I.3.1. Численное дифференцирование

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

, (11а)

, (11б)

, (11в)

Здесь – концентрация i-го реагента в момент времени t, рассчитанная при изменении константы скорости k_j на величину h, когда все остальные константы скорости сохраняют обычные значения. Главное достоинство такого метода — простота, поскольку здесь можно обойтись программой для решения прямой кинетической задачи.

Формулы (11) различаются по точности и требуют разного объема вычислительной работы. Так, в случае формулы 1-го порядка точности (11а) достаточно M раз решить прямую задачу, чтобы получить набор NM коэффициентов чувствительности. Более точная формула 2-го порядка (11б) требует 2M раз решить прямую задачу, а формула 4-го порядка (11в) — 4M раз.

Серьезный недостаток численного дифференцирования проявляется при оценке малых коэффициентов чувствительности. Действительно, если коэффициент мал, то даже значительное приращение k_j вызовет лишь малое изменение y_i. Поскольку точность решения прямой задачи ограничена заданной величиной  (см. п. I.2), то разность может оказаться меньше ошибки расчета, т. е. в лучшем случае удастся получить лишь грубую оценку порядка величины , а в худшем — результат будет совершенно бессмысленным.

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

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

I.3.2. Уравнения чувствительности

Дифференцируя кинетические уравнения (1) по константам скорости, можно получить систему дифференциальных уравнений для определения коэффициентов чувствительности:

.

Меняя порядок дифференцирования в левой части, окончательно получаем:

, (12)

где – элементы якобиана (5),
.

Очевидно, начальные условия для системы уравнений (12) имеют вид

,

поскольку начальные значения концентраций не зависят от констант скорости.

Уравнения (12) называются уравнениями чувствительности [1, 3]. Так как для вычисления элементов и в правой части необходимо знать концентрации реагентов в момент времени t, то уравнения (12) следует объединить с кинетическими уравнениями (1) в единую систему из N(M) дифференциальных уравнений. Численное интегрирование такой системы позволяет одновременно получать кинетические кривые и коэффициенты чувствительности. Этот подход называется прямым методом (ПМ) расчета чувствительностей [4].

ПМ часто применялся на практике, особенно в первых работах по анализу чувствительности для задач химической кинетики [1, 4]. Он обеспечивает достаточную степень точности независимо от величины коэффициентов . Основным недостатком ПМ является его трудоемкость: требуется решать систему дифференциальных уравнений большой размерности. Например, для 10‑стадийной модели высокотемпературного окисления азота с участием N₂, O₂, N, O и NO (см. п. I.4.2) необходимо интегрировать 55 дифференциальных уравнений.

I.3.3. Метод функций Грина

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