Пред.
След.
Макеты страниц
Распознанный текст, спецсимволы и формулы могут содержать ошибки, поэтому с корректным вариантом рекомендуем ознакомиться на отсканированных изображениях учебника выше Также, советуем воспользоваться поиском по сайту, мы уверены, что вы сможете найти больше информации по нужной Вам тематике ДЛЯ СТУДЕНТОВ И ШКОЛЬНИКОВ ЕСТЬ
ZADANIA.TO
В предыдущем параграфе было показано, как находить стационарные решения системы (5.1.1) при фиксированных значениях параметров. Если мы хотим изучить поведение динамической модели в целом, то обычно оказывается недостаточно знать ее характеристики только при одном конкретном значении того или иного параметра — нужно иметь представление о характере поведения модели в зависимости от значений параметров, изменяющихся в некотором диапазоне. В этом параграфе мы займемся построением зависимости стационарных решений от одного (скалярного) параметра $\alpha$, входящего в систему (5.1.3). Найденные зависимости, представленные в виде соответствующих графиков, мы будем называть диаграммой решений (см. гл. $3, \$ 3.1$ ). Простейший метод построения диаграммы заключается в последовательном использовании итерационной процедуры, описанной в §5.1. Так, например, выберем в интервале изменения значений параметра $\alpha$ узловые точки $\alpha^{0}<\alpha^{1}<\ldots<\alpha^{k}$ и применим последовательно метод Ньютона для указанных значений $\alpha$. В качестве начального приближения для нашего итерационного процесса при $\alpha=\alpha^{i}$ мы будем выбирать решение $\mathbf{x}$, найденное при $\alpha=\alpha^{i-1}$, т. е. $\mathbf{x}\left(\alpha^{i-1}\right)$. Если сетка значений параметра $\alpha$ достаточно густа и на полученных зависимостях отсутствуют точки бифуркации, то для обеспечения сходимости метода Ньютона при каждом значении $\alpha$ оказывается достаточным одной-двух итераций. Указанная методика не позволяет, однако, переходить через точки поворота на диаграмме решений (см. § 3.1), а в окрестности точек ветвления процесс может даже расходиться. Поскольку метод последовательных шагов не является универсальным, были развиты методики, с помощью которых построение зависимости решения от параметра в целом происходит более или менее автоматически. Основные принципы этих методик будут рассмотрены в последующих пунктах. 5.2.1. Отображение параметра Во многих случаях оказывается возможным использовать некоторые специальные свойства системы (5.1.3), которые позволяют разработать тот или иной неитерационный алгоритм (или же, при необходимости, итерационный, но в пространстве существенно меньшей размерности). Рассмотрим, например, случай, когда система (5.1.3) нелинейна лишь по одной переменной $x_{k}$ и линейна по всем другим переменным, а также по параметру $\alpha$. Тогда для построения диаграммы решений мы можем использовать следующий подход. Выберем последовательность значений переменной $x_{k}$ (достаточно близких друг к другу). Для каждого выбранного значения переменной $x_{k}$ решим систему линейных алгебраических уравнений (5.1.3) относительно неизвестных $x_{1}, x_{2}, \ldots, x_{k-1}, x_{k+1}, \ldots, x_{n}, \alpha$ (например, методом исключения Гаусса) — в результате мы получим одну из точек диаграммы решений. Существуют и другие возможности отображения параметра, когда, например, мы пользуемся тем, что умеем решать квадратное уравнение, можем построить соответствующую обратную функцию и т. д. Используемый при этом подход всегда зависит от конкретного вида уравнений, и мы продемонстрируем указанную методику, на нескольких примерах. Рассмотрим задачу 1 (см. гл. 4, п. 4.2.1). Стационарное состояние в данном случае описывается уравнениями где $x$ и $\boldsymbol{\Theta}$-переменные состояния. Умножая уравнение (5.2.1) на параметр $B$ и вычитая полученный результат из уравнения (5.2.2), после простых преобразований мы получаем соотношение которое позволяет свести систему двух уравнений (5.2.1) и (5.2.2) к одному нелинейному уравнению для переменной $\Theta$ : Уравнение (5.2.4) зависит от переменной $\Theta$ нелинейно, тогда как параметры $\mathrm{Da}^{1)}, B, \boldsymbol{\beta}$ и $\Theta_{c}$ входят в него линейным образом. Параметр $\Lambda$ (после умножения обеих частей уравнения на $\Lambda$ ) входит в полученное уравнение квадратичным образом, а параметр $\gamma$ — нелинейно. Для нахождения зависимости решения от параметра Da (при фиксированных значениях остальных параметров) можно использовать следующую процедуру: Пример одного из вариантов расчета представлен в табл. 5.3. Читатель может легко построить соответствующие алгоритмы для отображения параметров $B, \beta$ и $\Theta_{c}$, а если воспользоваться решением соответствующего квадратного уравнения, то и для отображения параметра $\Lambda$. Рассмотрим теперь задачу 1 в форме (P1-6), (P1-7), введя параметры $\tau$ (время задержки в реакторе) и а по формулам $\beta=a \tau, \mathrm{Da}=k_{0} \tau$. Совершенно аналогично можно построить ал- Таблица 5.3. Отображение параметра Dа в задаче $1(\gamma=20, B=10$, $\left.\Lambda=1, \Theta_{c}=-5, \beta=0,6\right)$. 1) горитм для отображения параметра $\tau$. Фиксируя значение $\Theta$, мы получаем для т квадратное уравнение вида Диаграмма решений, построенная с помощью этого уравнения для нескольких различных значений параметра $a$, представлена на рис. 5.2. Читатель может легко вывести аналогичные алгоритмы для задач 5 и 6. В обоих случаях фиксируются значения концентрации субстрата $C_{S}$ и подсчитывается значение некоторого параметра задачи. В задаче 6 это могут быть параметры $\hat{\mu}, V / F$, $K_{\text {s, }} K_{\text {i. }}$ В случае же задачи 5 соответствующая процедура оказывается более сложной, и мы приведем лишь краткое ее описание. 0. Левые части соотношений (P5-1)-(P5-6) полагаем равными нулю. При этом для стационарных значений $c_{z}$ и $c_{\text {т }}$ мы получим $c_{\mathrm{z}}=c_{\mathrm{z} 0}, c_{\mathrm{r}}=c_{\mathrm{т} 0}$. Рис. 5.2. Диаграмма стационарных решений задачи $1, \gamma=20, B=10, \Lambda=1$, $\boldsymbol{\theta}_{c}=-5, k_{0}=1$; сплошные линии — устойчивые решения, штриховые — неустойчивые решения. 4. Из уравнения (Р5-9) для $\mu$ (величину $\mu$ мы уже определили на этапе 2) можно найти теперь один из параметров $\hat{\mu}, K_{\mathrm{s}}, K_{\mathrm{i}}, K_{\mathrm{a}}, K_{1}$ в предположении, что остальные параметры заданы. Как уже говорилось, метод отображения параметра не обладает достаточной общностью. В связи с этим были развиты методы общего характера, которые мы рассмотрим в следующих пунктах. 5.2.2. Метод дифференцирования по параметру Пусть в точке ( $\mathbf{x}^{0}, \alpha^{0}$ ) выполнены условия теоремы о неявной функции: $\mathbf{f}\left(\mathbf{x}^{0}, \alpha^{0}\right)=\mathbf{0}$ и матрица Якоби $\mathbf{J}(\mathbf{x}, \alpha)=\left[\partial f_{i} / \partial x_{j}\right]$ невырожденная. Соотношение (5.2.6) можно трактовать как систему дифференциальных уравнений для нахождения функции $\mathbf{x}(\alpha)$ с начальным условием Если матрица $\mathbf{J}$ на промежутке $\left[\alpha^{0}, \alpha^{1}\right]$ оказывается регулярной (имеющей обратную), то зависимость решения от параметра $\mathbf{x}(\alpha)$, найденная путем интегрирования этой системы, будет удовлетворять соотношению поскольку В критических точках диаграммы решений матрица $\mathbf{J}$ оказывается вырожденной, и потому указанный подход, т. е. численное интегрирование системы дифференциальных уравнений (5.2.6), отказывает в тех же точках, что и процедура последовательного применения метода Ньютона 1). Однако этот метод можно модифицировать путем введения нового параметра. Обычно в качестве такого параметра выбирается длина дуги на кривой $\mathbf{f}(\mathbf{x}, \alpha)=\mathbf{0}$ в $(n+1)$-мерном пространстве переменных $x_{1}, x_{2}, \ldots, x_{n}, \alpha^{2)}$ Обозначим ее через $z$; тогда, дифференцируя тождество $\mathfrak{f}(\mathrm{x}(z), \alpha(z))=0$ по $z$, находим При этом уравнение определяет параметр $z$ как длину дуги кривой. Начальное условие (5.2.10) теперь принимает вид Соотношения (5.2.8) можно рассматривать как систему $n$ линейных алгебраических уравнений относительно $n+1$ неизвестных $d x_{1} / d z, d x_{2} / d z, \ldots, d x_{n} / d z, d \alpha / d z$. Мы можем решить эту систему так, чтобы найти зависимость выбранных $n$ неизвестных от заданной неизвестной $d x_{k} / d z$, где $k$ фиксировано. В дальнейшем для упрощения записи будем обозначать Для построения указанной зависимости необходимо, чтобы матрица которая получается из матрицы системы (5.2.8) $\tilde{\mathbf{J}}$ вычеркиванием $k$-го столбца (квадратная матрица размером $n \times n$ ), оказалась невырожденной. В точке поворота всегда можно выбрать индекс $k$ так; чтобы матрица $\mathbf{J}_{k}$ была невырожденной; в точке ветвления, наоборот, все $\mathbf{J}_{k}$ вырождены. Если матрица $\mathbf{J}_{k}$ — невырожденная, то систему (5.2.8) можно представить в виде где коэффициенты $\beta_{i}$ подсчитываются методом исключения Гаусса. Выбор индекса $k$ представляет собой самостоятельную проблему. Мы либо фиксируем его заранее, либо видоизменяем алгоритм исключения Гаусса с выбором главного элемента для решения (неквадратной) системы $n$ линейных алгебраических уравнений (5.2.12) относительно $n+1$ неизвестных. При этом неизвестная, в столбце которой отсутствует главный элемент, есть неизвестная $x_{k}{ }^{11}$. Подставим найденные зависимости (5.2.13) в уравнение (5.2.9) : Таким образом, мы определили производную $d x_{k} / d z$ с точностью до знака. Если индекс $k$ остается неизменным для всего продолжения, то знак в формуле (5.2.14) в процессе продолжения не меняется 1). Отсюда следует, что тогда $x_{k}$ во время всего процесса продолжения должно либо постоянно возрастать, либо постоянно убывать. Однако это не всегда бывает так, и тогда фиксированный выбор $k$ невозможен. Поэтому в практических ситуациях выбор $k$ осуществляется адаптивным образом — значение индекса $k$ получается в результате применения метода исключения Гаусса (см. выше). При этом знак в формуле (5.2.14) будет определяться тем, возрастает или убывает в этот момент выбранная величина $x_{k}$ вдоль кривой в направлении возрастания $z$. Поэтому мы введем «знаковые» переменные $N_{i}$, положив Указанные «параметры направления» будут вычисляться на каждом шаге процесса продолжения. Выбор знака в формуле (5.2.14) определяется, следовательно, соответствующей величиной $N_{k}$, взятой с предыдущего шага. Отметим, что оставшиеся производные $d x_{i} / d z \quad(i Проиллюстрируем использование описанного метода продолжения на примере нахождения зависимости стационарного решения задачи 1 от времени задержки $\tau$ (ср. рис. 5.2). Уравнения (5.1.3) в указанном случае переписываются в виде Начальное условие (5.2.10) при $\tau=0$ принимает вид Рис. 5.3. Схема метода дифференцирования по длине дуги и метода продолжения Эйлера. 1) При выборе «главного» элемента в строке 1 берется $\max _{j}\left(p ; a_{1 j}\right)$. Аналогично на следуюших шагах метода исключения. a (расширенная) матрица Якоби $\widetilde{\mathbf{J}}$ системы (5.2.16) (если обозначить $E=\exp (\Theta /(1+\Theta / \gamma)))$ записывается как Применим теперь вычислительный алгоритм, соответствующий расчетной схеме рис. 5.3, для случая различных шагов интегрирования $\Delta z$. Результаты вычислений представлены в табл. 5.4. Первая часть этой таблицы может быть использована для более глубокого уяснения метода. Читатель может легко сопоставить полученные результаты с данными рис. 5.2 для $a=1$. Вторая часть таблицы характеризует влияние погрешностей аппроксимации при выбранном методе интегрирования. При этом чем меньше шаг интегрирования, тем лучше найденные значения аппроксимируют реальную зависимость. Из таблицы видно, что даже в случае рассматриваемой простой задачи для получения достаточно точных результатов нам пришлось бы выбирать $\Delta z$ очень малым. Это обусловлено, конечно, в первую очередь низким порядком аппроксимации используемого метода Эйлера (погрешность аппроксимации в этом случае есть $O(\Delta z)$ ). Если выбрать более эффективный метод интегрирования, например метод Рунге — Кутты 4-го порядка, то результаты окажутся более точными. Полученное решение, однако, также, хотя и в меньшей степени, будет отличаться от истинной зависимости. Вследствие этого представляется необходимым после прохождения некоторого интервала изменения переменной $z$ всегда уточнять решение, т. е. корректировать его так, чтобы выполнялось исходное уравнение (5.1.3), или, для рассматриваемого случая, уравнения (5.2.16). Так появились методы продолжения типа предиктор-корректор, которые мы опишем в следующем пункте. 5.2.3. Алгоритм продолжения типа предиктор-корректор В принципе существуют две возможности коррекции погрешностей аппроксимации, накапливающихся при продолжении решения с использованием простого предиктора. Первая из них заключается в том, чтобы контролировать, как отличается решение от истинного, например, по норме $\|\mathbf{f}\|$, и когда это отличие превышает некоторый заранее заданный допуск, использовать для уточнения результатов какой-либо итерацион- ный метод (например, метод Ньютона). При другом подходе также используется метод Ньютона, но после каждого шага предиктора. При этом в методе Ньютона осуществляется лишь такое число итераций, чтобы приращение $\Delta \mathbf{x}$ по соответствующей норме оказалось бы меньше заданной точности вычислений. В обоих подходах уточняются выбранные $n$ из $n+1$ переменных $x_{1}, \ldots, x_{n}, x_{n+1}=\alpha$. Выбор неизменяемой переменной (ее индекса $k$ ) не является произвольным: например, в окрестности предельной точки это не может быть $\alpha$, а в окрестности экстремума $x_{i}$ это не может быть $x_{i}$. Поэтому алгоритм выбирает ее сам, причем одной из возможностей является использование того же самого механизма, что в п. 5.2.2. Диаграмма вычислений для описанного алгоритма продолжения типа предиктор-корректор представлена на рис. 5.4. В работе [5.5] этот алгоритм расписан на языке Фортран, соответствующая подпрограмма называется DERPAR. В данной книге мы будем использовать это название. Отметим, что в этой подпрограмме вместо метода Эйлера в предикторе используются варианты многошагового метода АдамсаБашфорта переменного порядка. Приведенный алгоритм продолжения беспрепятственно преодолевает точки поворота на зависимостях решения от параметра. В точках же ветвления он обычно дает продолжение первоначальной ветви решения. Продолжение может не удаться, если очередная вычисленная точка окажется слишком близкой к точке ветвления. Сигналом, указывающим, что мы перешли за точку ветвления, может служить изменение знака $\operatorname{det} \mathbf{J}_{k}$ определителя матрицы $\mathbf{J}_{k}$. Таким образом, указанный алгоритм формирует зависимость решения от параметра, представляющую собой непрерывную кривую в $(n+1)$-мерном пространстве переменных $x_{1}, \ldots, x_{n}, x_{n+1}=\alpha$. Разумеется, алгоритм DERPAR не приспособлен для того, чтобы с его помощью построить всю диаграмму решений целиком. Для каждой изолированной ветви кривой $f(x, \alpha)=0$ ему нужна новая начальная точка, которую можно получить, например, методом случайного перебора начальных приближений для схемы Ньютона (см. § 5.1). Наконец, отметим, что описанный алгоритм не рассчитан на точное нахождение бифуркационных точек в пространстве $(x, \alpha)$. Обнаружив такую точку, ее координаты можно точно вычислить с помощью алгоритмов, описанных в п. 5.4.1, и определить затем соответствующие начальные значения для процесса продолжения на остальных ветвях решения (см. п. 5.4.2). В последнее время был разработан целый ряд методов и алгоритмов, связанных с процессом продолжения (см., например, работы [5.24], [5.25]), где большей частью используется практически тот же подход, что и в алгоритме DERPAR. Результат применения алгоритма продолжения DERPAR для задачи 2 [5.5] представлен на рис. 5.5. При этом вся кривая целиком была построена путем одного обращения к алгоритму. Отметим, что на диаграмме решений имеется 6 точек поворота и отсутствуют точки ветвления. Точка пересечения трех ветвей на рисунке не является точкой бифуркации, по- Рис. 5.5. Диаграмма стационарных решений задачи $2, \gamma=1000, B=22$, $\boldsymbol{\beta}_{1}=\boldsymbol{\beta}=2, \Theta_{c_{1}}=\Theta_{c_{2}}=0: \boldsymbol{\Lambda}=1, \boldsymbol{\alpha} \equiv \mathrm{Da}_{1}=\mathrm{Da}_{2} ;$ сплошные линии — устойчивые решения, штриховые — неустойчивые решения. Рис. 5.6. Диаграммы стационарных решений некоторых задач из главы 4; сплошные линии — устойчивые решения, штриховые линии — неустойчивые решения. a) Задача 8, каскад из двух реакторов с однонаправленным течением, уравнения (P8-2a), (P8-3a), (P8-4), (P8-5); $A=2, B=4, D_{2}=$ $=10 D_{1}$. b) Задача 4; $\gamma=3, \beta=1,5 ; v_{0}=0,01$. c) Задача $10 ; \sigma=16$, $b=4 . d$ ) Задача 8; $N=2$, уравнения (P8-2)-(P8-5); $A=2, B=6, D_{2}=$ $\left.=10 D_{1} . e\right)$ Задача 8; $N=3$, уравнения (P8-7)-(P8-12); $A=2, B=6$. $D_{2}=10 D_{1} . f$ ) Задача 8; $N=4$, уравнения (P8-7)-(P8-10), (P8-14)-(P8-17), $\left.A=2, B=6, D_{2}=10 D_{1} . g\right)$ Задача $2 ; \gamma=20, \Lambda=1, \Theta_{c_{1}}=\theta_{c_{2}}=-5$, $\beta_{1}=\beta_{2}=1, B=15, \mathrm{Da}=\mathrm{Da}_{1}=\mathrm{Da}_{2}$. h) задача $1 ; \gamma=20 ; B=9, \boldsymbol{\theta}_{c}=0$, $\mathrm{Da}=0,02$. $i$ ) Задача 5 ; значения параметров см. в табл. 4.1. Построение диаграммы решений представляет собой важнейшую задачу анализа нелинейной динамической системы1). Читателю, который хочет практически освоить эту проблематику, мы советуем рассчитать несколько диаграмм решений самостоятельно. На рис. 5.6 изображены 9 диаграмм решений для различных задач, рассмотренных в гл. 4. Диаграммы решений, представленные на рис. $5.6 \mathrm{~b}, \mathrm{~h}, \mathrm{i}$, можно получить с помощью отображения соответствующего параметра. Диаграмму на рис. 5.6c можно построить аналитически. Начальные значения для продолжения решений на рис. $5.6 \mathrm{~d}$, e, f читатель найдет в табл. 5.1 при $D_{1}=0,4$.
|
1 |
Оглавление
|