d -ш

7<: б<; 5<:

з<; 2

о

- i I- -1-

10 20 30- ;С 50 60 70 80 90 100 ПО 120 130 М

К

Рис. 5.1. Диаграмма успешности применения метода Ньютона в зависимости' от случайно выбранного начального приближения. К-порядковый номер-начального приближения. На вертикальной оси приведены номера полученных решений (см. табл. 5.1 при = 4). Символ d означает расходимость метода от соответствующего начального приближения. Для решений, отмеченных в табл. 5.1 звездочкой, два симметричных решения на рисунке представлены

отдельно.

Проиллюстрируем указанный подход на примере, приведенном в гл. 4 в качестве задачи 8 для случая = 2,3,4. В табл. 5.1 представлены решения для двух, трех и четырех связанных между собой реакторов (т. е. для п = 4, 6, 8 в системе (5.1.3)). Число решений при этом оказывается соответственно равным 3, 7 и 15 [5.2]. В случае цепочки из большего числа реакторов количество различных стационарных решений будет еще больше.

Начальные приближения х° для метода Ньютона выбирались случайным образом из интервалов

хОе(0, 10), г = 1, 2, .. ., п. (5.1.8)

Примеры расчета итераций приведены в табл. 5.2. На рис. 5.1 изображена диаграмма успешности этого подхода для случая четырех связанных между собой реакторов {N = 4).

выбранного начального приближения мы практически получаем один из следующих результатов:

1. итерационный процесс сходится к физически приемлемому решению (новому или уже найденному ранее);

2. итерационный процесс сходится к физически неприемлемому решению;

3. итерационный процесс расходится, т. е. значения переменных с какого-то момента в ходе итераций превышают установленные пределы;

4. итерационный процесс не сходится с предписанной точностью при числе итераций, не превосходящем заданного.

Таблица 5.1. Стационарные решения в случае модели линейного каскада связанных друг с другом реакторов {задача 8: Л = 2, S = 6, Di = 0,4, Dz = 4). В первой строке каждой ячейки указана величина концентрации X, во второй - величина коицентраЩ1и Y.

Каждое решение, отмеченное звездочкой, задает по существу два решения. Второе решение, симметричное к указанному в таблице, например для N=4, мы получаем, меняя местами значения концеитраци* в первом н в четвертом, а также во втором н в третьем реакторах.

Таблица 5.2. Поведение итераций в процессе вычисления стационарных решений методом Ньютона (задача 8: = 4, Л = 2, S = 6, Z), = 0,4, Ог = 4).

Область сходимости, т. е. область подходящих начальных приближений, в некоторых случаях может оказаться очень мала. Иногда в подобной ситуации используется метод продолжения решения из области (характеризуемой определенным значением выбранного параметра), в которой решение отыскать сравнительно легко, в область, где нахождение решения представляет определенные трудности. Параметр, по которому мы продолжаем решение, может быть реальным (физическим) параметром задачи, либо может искусственно вводиться в постановку задачи. (Иногда мы говорим о погружении исходной задачи в класс задач, зависящих от параметра.) Первый из этих подходов будет рассмотрен в § 5.2, а второй, так называемый метод введения искусственного параметра, кратко описывается ниже.

Запишем систему (5.1.3) в векторной форме

f(x) = 0; (5.1.9)

при этом значение параметра а мы считаем фиксированным и поэтому в уравнениях его не указываем. Идея метода введения искусственного параметра заключается в построении функции Н(т, х), зависящей от некоторого вещественного параметра х, причем обычно те [0,1]. Функция Н выбирается таким образом, чтобы решение уравнения

Н(0, х) = 0 (5.1.10)

можно было построить по возможности просто, например, не прибегая к процессу итераций (это решение мы обозначим через Хо). Далее, уравнение

Н(1, х) = 0 (5.1.11)

ДОЛЖНО иметь те же корни, что и уравнение (5.1.9), например

Н(1, x) = f(x). (5.1.12>

Если существует непрерывная зависимость х(т), описывающая решение уравнения

Н(т, х(т)) = 0, те [0,1], (5.1.13)

с условием х(0) = хо, то х(1) является решением исходного уравнения (5.1.9). В принципе мы можем протянуть решение X от т = 0 (когда х = Хо) до т = 1 двумя способами. Первый из НИХ заключается в последовательном использовании какого-нибудь итерационного метода (например, метода Ньютона) в узлах сетки то = О, ti, тг, Тт = 1. В качестве начального приближения процесса итераций для т = т/ можно выбрать решение x(t/ i). Второй способ основан на дифференцировании тождества (5.1.13) по т, т. е. зависимость х(т) ищется как решение системы дифференциальных уравнений

- = -[Н;(г. х)]-н;(г, X) (5.1.14)

с начальным условием х(0) = Хо.

Укажем две возможности выбрать функцию Н:

Н(т, x) = f(x)-b(T-l)f(xo) (5.1.15)

Н(т, x) = (l-T)(x-Xo)-bTf(x). (5.1.16)

Представление (5.1.15) тесно связано с методом Ньютона. (Чтобы убедиться в этом, нужно, взяв Н в виде (5.1.15), решить дифференциальное уравнение (5.1.14) методом Эйлера с шагом h=l.) Более подробный анализ различных вариантов метода введения искусственного параметра читатель может найти в. книгах [5.1, 5.3]. Описание некоторых других подходов и их связь с итерационными методами рассмотрены в работе [5.4].

5.2. ЗАВИСИМОСТЬ СТАЦИОНАРНЫХ РЕШЕНИИ ОТ ПАРАМЕТРА-ДИАГРАММА РЕШЕНИЙ

В предыдущем параграфе было показано, как находить стационарные решения системы (5.1.1) при фиксированных значениях параметров. Если мы хотим изучить поведение динамической модели в целом, то обычно оказывается недостаточно знать ее характеристики только при одном конкретном значении того или иного параметра -нужно иметь представление о характере поведения модели в зависимости от значений параметров, изменяющихся в некотором диапазоне.

В ЭТОМ параграфе мы займемся построением зависимости стационарных решений от одного (скалярного) параметра а, входящего в систему (5.1.3). Найденные зависимости, представленные в виде соответствующих графиков, мы будем называть диаграммой решений (см. гл. 3, § 3.1).

Простейший метод построения диаграммы заключается в последовательном использовании итерационной процедуры, описанной в § 5.1. Так, например, выберем в интервале изменения значений параметра а узловые точки а <; а'< ... <; а* и применим последовательно метод Ньютона для указанных значений а. В качестве начального приближения для нашего итерационного процесса при а = а' мы будем выбирать решение х, найденное при а = а'~, т. е. х(а'-). Если сетка значений параметра а достаточно густа и на полученных зависимостях отсутствуют точки бифуркации, то для обеспечения сходимости метода Ньютона при каждом значении а оказывается достаточным од-ной-двух итераций. Указанная методика не позволяет, однако, переходить через точки поворота на диаграм.ме решений (см. § 3.1), а в окрестности точек ветвления процесс может даже расходиться. Поскольку метод последовательных шагов не является универсальным, были развиты методики, с помощью которых построение зависимости решения от параметра в целом происходит более или менее автоматически. Основные принципы этих методик будут рассмотрены в последующих пунктах.

5.2.1. Отображение параметра

Во многих случаях оказывается возможным использовать некоторые специальные свойства системы (5.1.3), которые позволяют разработать тот или иной неитерационный алгоритм (или же, при необходимости, итерационный, но в пространстве существенно меньшей размерности). Рассмотрим, например,случай, когда система (5.1.3) нелинейна лишь по одной переменной Xk и линейна по всем другим переменным, а также по параметру а. Тогда для построения диаграммы решений мы можем использовать следующий подход. Выберем последовательность значений переменной Xk (достаточно близких друг к другу). Для каждого выбранного значения переменной Xk решим систему линейных алгебраических уравнений (5.1.3) относительно неизвестных лгь лгг, х^-и Xk+u , Хп,а (например, методом исключения Гаусса) - в результате мы получим одну из точек диаграммы решений. Существуют и другие возможности отображения параметра, когда, например, мы пользуемся тем, что умеем решать квадратное уравнение, можем построить соответствующую обратную функцию и т. д. Используемый при этом

Напомним, что Da - обозначение одного параметра (числа Дамкёлера), а не произведение D-a. - Прим. ред.

ПОДХОД всегда зависит от конкретного вида уравнений, и мы продемонстрируем указанную методику на нескольких примерах.

Рассмотрим задачу 1 (см. гл. 4, п. 4.2.1). Стационарное состояние в данном случае описывается уравнениями

-Ax-bDa(l-x)exp(jq) = 0, (5.2.1)

-Ав + ВВг{1-х) ехр (iq) - Р(в - в,) = О, (5.2.2)

где X и в - переменные состояния. Умножая уравнение (5.2.1) на параметр В и вычитая полученный результат из уравнения (5.2.2), после простых преобразований мы получаем соотношение

= Т + (5.2.3)

которое позволяет свести систему двух уравнений (5.2.1) н (5.2.2) к одному нелинейному уравнению для переменной в:

-Ав-Ь Da (в - в --f) ехр (- Р (0-Ь 6,) = 0.

(5.2.4)

Уравнение (5.2.4) зависит от переменной в нелинейно, тогда как параметры Da В, р и вс входят в него линейным образом. Параметр Л (после умножения обеих частей уравнения на Л) входит в полученное уравнение квадратичным образом, а параметр у - нелинейно. Для нахождения зависимости решения от параметра Da (при фиксированных значениях остальных параметров) можно использовать следующую процедуру:

1. Выбираем значение в(в > вс)..

2. Из уравнения (5.2.4) подсчитываем значение параметра Da.

3. Из уравнения (5.2.3) находим значение другой переменной состояния X.

Пример одного из вариантов расчета представлен в табл. 5.3. Читатель может легко построить соответствующие алгоритмы для отображения параметров В, Р и вс, а если воспользоваться решением соответствующего квадратного уравнения, то и для отображения параметра Л.

Рассмотрим теперь задачу 1 в форме (Р1-6), (Р1-7), введя параметры х (время задержки в реакторе) и а по формулам Р = ах, Da = oT. Совершенно аналогично можно построить ал-

Таблица 5.3. Отображение параметра Da в задаче 1 {у = 20, В Л= 1, е. = -5, Р = 0,6).)

= 10,

1) Cm. п. 5.3.4 (S - седло, SU - устойчивый узел, NU - неустойчивый узел, SO - устойчивый фокус, NO - неустойчивый фокус. Я ,s - собственные числа матрицы Якоби).

горитм ДЛЯ отображения параметра т. Фиксируя значение в, мы получаем для т квадратное уравнение вида

- kpa (е - вс) -- / 0 Л

-expl

(тТвл)] +

+ [к, (В - в) ехр (-npes) - ° (® - 9.)] + [-Лв] - 0. (5.2.5)

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

Читатель может легко вывести аналогичные алгоритмы для задач 5 и 6. В обоих случаях фиксируются значения концентрации субстрата Cs и подсчитывается значение некоторого параметра задачи. В задаче 6 это могут быть параметры р V/F, Ks, Ki. В случае же задачи 5 соответствующая процедура оказывается более сложной, и мы приведем лишь краткое ее описание.

0. Левые части соотношений (р5-1) -(р5-6) полагаем равными нулю. При этом для стационарных значений Cz и Ст

мы получим Cz = CzO, Ст = Сто-

1. Выбираем значение Cs.

2. Из уравнений (р5-1) и (р5-2) вычисляем Сх и р.

3. Подставляя результат в зфавнение (р5-3), находим линейную зависимость рс от Сс, которую затем подставляем

Рис. 5.2. Диаграмма стационарных решений задачи 1, у = 20, В = 10, Л=1, вс = -5, = 1; сплошные линии - устойчивые решения, штриховые - неустойчивые решения.

в соотношение (р5-4). При этом мы получаем квадратное уравнение относительно Сс. Физически допустимые решения этого квадратного уравнения дают нам окончательные значения Сс, рс-

4. Из уравнения (р5-9) для р, (величину р, мы уже определили на этапе 2) можно найти теперь один из параметров р.,/Cs, Ки /Са, К\ в предположении, что остальные параметры заданы.

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

При этом уравнение

т+-+т+т= < ->

Точнее, в этих точках теряет смысл само уравнение (5.2.6). -Прим.

> Дело, конечно, не в выборе параметра, а в рассмотрении х и а как равноправных переменных. Эта, почти очевидная, мысль не требует обязательного рассмотрения дифференциального уравнения (5.2.6). - Прим. ред.

5.2.2. Метод дифференцирования по параметру

Пусть в точке (х°, а ) выполнены условия теоремы о неявной функции: f(x°,a°) = 0 и матрица Якоби J{x, a) = [dfi/dxj] невырожденная.

Тогда для функции х{а), задаваемой уравнением f(x, а) = 0, g = -[J(x, а)Г'{х, а). (5.2.6)

Соотношение (5.2.6) можно трактовать как систему дифференциальных уравнений для нахождения функции х(а) с начальным условием

х(ао) = х<. (5.2.7)

Если матрица J на промежутке [а°, а'] оказывается регулярной (имеющей обратную), то зависимость решения от параметра х(а), найденная путем интегрирования этой системы, будет удовлетворять соотношению

f(x(a), а) = 0, аеК, а'],

поскольку

-£ f (X (а), а) = J (х (а), ) + (х (а), а) = 0.

В критических точках диаграммы решений матрица J оказывается вырожденной, и потому указанный подход, т. е. численное интегрирование системы дифференциальных уравнений (5.2.6), отказывает в тех же точках, что и процедура последовательного применения метода Ньютона Однако этот метод можно модифицировать путем введения нового параметра. Обычно в качестве такого параметра выбирается длина дуги на кривой f(x,a) = 0 в (п-f-1)-мерном пространстве переменных

Xi,X2.....л: , Обозначим ее через г; тогда, дифференцируя

тождество f(x(z),a(z)) = 0 по г, находим

определяет параметр z как длину дуги кривой. Начальное условие (5.2.10) теперь принимает вид

2 = 0: х = х , а = аР. (5.2.10)

Соотношения (5.2.8) можно рассматривать как систему п линейных алгебраических уравнений относительно tt-f-1 неизвестных dxi/dz, dx2/dz.....dxn/dz, da/dz. Мы можем решить

эту систему так, чтобы найти зависимость выбранных п неизвестных от заданной неизвестной dx/dz, где k фиксировано. В дальнейшем для упрошения записи будем обозначать

(5.2.11)

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

Г dfi df, dfi df, -I

J* =

udx,

dfn dx

dfn dx

n+l -I

(5.2.12)

которая получается из матрицы системы (5.2.8) J вычеркиванием k-ro столбца (квадратная матрица размером пУп), оказалась невырожденной.

В точке поворота всегда можно выбрать индекс k так, чтобы матрица Jk была невырожденной; в точке ветвления, наоборот, все Jk вырождены. Если матрица ik - невырожденная, то систему (5.2.8) можно представить в виде

dx

i=\, 2,

-bl, (5.2.13)

где коэффициенты p, подсчитываются методом исключения Гаусса. Выбор индекса k представляет собой самостоятельную проблему. Мы либо фиксируем его заранее, либо видоизменяем алгоритм исключения Гаусса с выбором главного элемента для решения (неквадратной) системы п линейных алгебраических уравнений (5.2.12) отнооительно п-\-\ неизвестных. При этом неизвестная, в столбце которой отсутствует главный элемент, есть неизвестная Xk

Подставим найденные зависимости (5.2.13) в уравнение (5.2.9):

г п + 1 -,-1/2

1 + (5.2.14)

i = l

- JV-fe -1

> Речь идет о выборе главного элемента по строке . Прн этом переменные могут & ть предварительно масштабированы, см. [5.10]. - Лрцл. ред.