Точный смысл сказанного таков. Разрешив меняться Р, мы получаем систему fi(T]i, т]2, Da, Р) = О (/ = 1, 2, 3). Ее решение определяет кривую в пространстве (rii, TI2, Da, Р), связные компоненты которой можно находить используя какой-либо алгоритм продолжения . Проекция найденной кривой на плоскость (Da, Р) дает часть линий бифуркационной диаграммы. - Прим. ред.

Мы называем эти точки точками вторичной комплексной бифуркации. 21 М. Холодннок н др.

<6.2.32), (6.3.27), (6.3.28), так же как в п. 5.4.1, используется метод Гаусса - Ньютона. Приложение указанного подхода к задаче И читатель может найти в работе [5.10].

Продемонстрируем использование этого метода для нахождения точек поворота в задаче 14 в случае, когда параметром является число Da. Таким образом, нас будет интересовать нахождение y{z), 6(2)-решения краевой задачи (6.1.35), (6.1.36). Недостающие начальные условия в точке z = 1 выберем как в (6.1.37): г/(1) = т11, 6(1) = TI2. Требуя выполнения краевых условий в точке z = 0 (см. 6.1.36а), мы получаем два уравнения:

РЛПи Ti2. Da) Л Ремг/(0)-/(0) = 0, Fim, TI2, Da)PeH0(O)-e(O) = O

(см. 6.2.41). Третье, бифуркационное уравнение здесь имеет вид (ср. с (6.2.42)):

F, (г,. Da) = (Рер^, (0) - р^, (0)) (Pep, (О - р^ (0)) -

- (РбнРв! (0) - Рш (0)) (РемР.2 (0) - Р;2 (0)) = О- (6-3.29)

Результаты, полученные при решении этих 3 уравнений методом Ньютона (с использованием уравнений второго порядка в вариациях), приведены в табл. 6.10. Читатель может сравнить координаты полученных точек поворота с данными рис. 6.8. Если продолжить теперь полученные координаты точек поворота в зависимости от какого-либо другого параметра задачи, то мы получим соответствующую бифуркационную диаграмму. Пример такого продолжения по параметру представлен на рис. 6.12 в плоскости параметров Da - р [6.18]

6.3.3. Нахождение точек комплексной бифуркации (бифуркации Андронова - Хопфа)

Нахождение точек комплексной бифуркации для тривиального пространственно-однородного решения рассматривалось нами в п. 6.3.1. В этом пункте мы займемся нахождением точек комплексной бифуркации на нетривиальных ветвях стационарных решений 2). Наиболее простым методом является здесь

Таблица 6.10. Метод Ньютона для определения точек поворота в задаче 14 (Y = 20, В = 15, р = 2, вс = О, Ре = 10, Рен = б)

метод имитации динамического поведения для последовательности значений исследуемого параметра. Начальные условия при этом выбираются в достаточно малой окрестности стационарного-решения. При таком подходе можно обнаружить потерю устойчивости стационарного решения в процессе продолжения этого решения по параметру. Ясно, что указанным способом можнО' находить только те точки комплексной бифуркации, в которых, (при продвижении вдоль ветви стационарных решений) стационарное решение теряет устойчивость.

S.3.3.I. Метод прямых - переход к системе с сосредоточенными параметрами

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

dt =тф(-1~2 + +) + ((6-3.31а) Граничные условия (4.3.9) дают

x,it)x, xjt)x, ying, yjt)g. (6.3.31b)

Таким образом, мы имеем систему 2{п-1) обыкновенных дифференциальных уравнений, для решения которой можно использовать методы, рассмотренные в § 5.5. Система (6.3.31а, Ь) представляет собой некоторую аппроксимацию исходной распределенной системы. Погрешность аппроксимации зависит от h и имеет порядок 0{h). При уменьшении h (т. е. с увеличением точности аппроксимации) возрастают размеры системы (6.3.31) и возникают трудности с вычислением собственных чисел соответствующих матриц. Поэтому описываемый подход применим лишь для не очень точного нахождения точек комплексной бифуркации, хотя в некоторых случаях даже при сравнительно большом шаге h ои дает вполне удовлетворительные результаты. Дискретизацию по координате z можно проводить и с помощью более точных формул, например посредством ортогональной коллокаций. Читателя, который более глубоко заинтересуется применением метода прямых для нахождения точек комплексной бифуркации, мы отсылаем к оригинальным работам [6.19, 6.20, 6.21, 6.22].

6.3.3.2. Прямой итерационный метод

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

альных уравнений, методы анализа которых рассматривались в гл. 5.

Рассмотрим, к примеру, систему типа реакция - диффузия (4.3.7) с граничными условиями типа 1 (см. (4.3.9)). Выберем равномерную сетку узловых точек (zt-ih, 1 = 0, 1, п, h = \/п) и обозначим приближенные значения решения в этих точках как

XI (t) - X (2 /), г/, (/) ~ у (2ь t). (6.3.30)

Аппроксимируем теперь производные по координате z с помощью простейших (трехточечных) разностных формул. Мы получим систему обыкновенных дифференциальных уравнений для t = 1, 2, п- 1:

dx D

зации) позволяет находить точку комплексной бифуркации гораздо точнее. Результаты, полученные с помощью первых двух подходов, удобно использовать в качестве начального приближения для прямого итерационного метода. Обратимся вновь к системе (4.3.7) с граничными условиями типа 1, т. е. условиями x{0,t)= x{\,t) = x\ y{0,t) = y{l,t)g (см. (4.3.9)). Обозначим через S линеаризованный оператор правых частей уравнений (4.3.7) и вычислим его для стационарного решения Xq{z), уо{г), т. е. для решения уравнений (6.1.1), (6.1.2). Вид

оператора задается выражением (6.3.2) Оц =(0(2), Уо()) и т.

В точке комплексной бифуркации оператор 2 имеет чисто мнимое собственное число К = is и, следовательно, имеется ненулевое решение и уравнения

Sa = isu; u

(6.3.32)

Положим u = Ur -f tUi, V = Vr + ivi И w = wr + iwi. Отделяя

теперь вещественную и мнимую части в соотношении (6.3.32), находим

Svlr = -SU (6.3.33) Sni = SUr,

или, более подробно.

, [df , df

VI +

= -svi.

(6.3.34)

I dg Vi + Wi

= -SWu

= SWr.

Граничные условия для функций vr, vi, Wr, Wi с учетом формул (4.3.9) имеют вид

Vr (0) = Vi (0) = Wr (0) = Wi (0) = 0, (6.3.35)

Vr (1) = Vi (1) = WR (1) = tWi (1) = 0. (6.3.36)

Собственная функция u = {v, w) определена с точностью до ненулевого комплексного множителя. Это означает, что для функций v{z) и w{z) мы можем достаточно произвольно задать два фиксированных ненулевых значения. Двумя возмож-

ными примерами такого выбора могут служить условия

1>к(1) = 1, ад(4-) = 1 (6.3.37а)

или

yjj(0)=l, w[{0)=l. (6.3.37b)

Читатель сам может предложить другие комбинации.

Таким образом, мы получили систему обыкновенных дифференциальных уравнений (6.1.1), (6.1.2), (6.3.34) (суммарно 12-го порядка) с 14-ю граничными условиями (4.3.9), (6.3.35) - (6.3.37) и с двумя параметрами s, L, которые нам предстоит определить в точке комплексной бифуркации. Число уравнений согласовано с числом неизвестных, и описанную задачу можно пытаться решать либо с помощью разностных методов, либо методом стрельбы. Рассмотрим здесь оба этих подхода.

Разностный подход. Заменим вторые производные простейшими трехточечными центрально-разностными соотношениями, В результате вместо уравнений (6.1.1), (6.1.2) и (6.3.34) мы получим следующую систему (ср. с (6.3.31)):

Xi i - 2л:, + л:,+, + = О, f= f (л:,-, у,),

121,2

yi-i-2yi + yi+i+-Ogi = 0, gi = g{Xi,yi), (6.3.38)

I2f2

f r. -l - 2t)R, , -f t)R, ,- + 1 -f (fx ,-t)R, i -f fy, iWR. i) = -SVl, i,

l 2fj2

wi, .--i - 2wi, i -f ayj. i+i + {gx, if i, i + gy. ii. i) = swr, i,

где t = 1, ..., П - 1.

Дополним эту систему очевидными граничными условиями:

Ха = Хп = X, уд = уп = у, Vr, о = ui, о = аУк. о = Wi, о = 0;

Ur. n=vi,n = wr,n = Wi, n = 0.

Тем самым мы получим систему 6(п--1) нелинейных (алгебраических) уравнений относительно 6{п-\-1)-{-2 неизвестных

Хи Уи fR, Vi,i, WR,iWi,i, (/ = 0, п), s, L.

В качестве двух недостающих условий можно, например, взять соотношения (6.3.37а). Полученную в результате систему 6(п--1)--2 нелинейных алгебраических уравнений можно решать с помощью метода Ньютона.

В п. 6.3.3.1 мы рассмотрели аппроксимацию исходных дифференциальных уравнений в частных производных системой

обыкновенных дифференциальных уравнений с помощью метода прямых. Точку комплексной бифуркации для этой аппроксимирующей системы (6.3.31) можно найти с помощью одного из методов, описанных в параграфе 5.5. При этом, если воспользоваться для решения указанной системы прямым итерационным методом в его несокращенной версии, которая была описана в п. 5.5.3, то нетрудно показать, что получающаяся нелинейная (алгебраическая) система имеет тот же самый вид, что и система (6.3.38), если подходящим образом выбрать для нее граничные условия Это утверждение читатель может легко проверить самостоятельно.

Метод стрельбы. Этот подход мы продемонстрируем на примере задачи 14, т. е. уравнений (Р14-7), (Р14-8) с граничными условиями (Р14-9), (Р14-10). Правые части этих уравнений отличаются от правых частей уравнений (4.3.7) членами с первой производной; кроме того, имеются различия и в граничных условиях (представляющих собой комбинацию граничных условий 1-го и 2-го рода). В качестве бифуркационного параметра выберем параметр Da и положим Le = 1.

Обозначим Е = {\ - i/)exp (в/(1-f в/)) и запишем уравнения, которым должны удовлетворять стационарные решения

1 .y y + Da£ = 0,

м

e -e + BDa£-p(e-e,) = 0.

(6.3.39)

Граничные условия при этом имеют вид

г/(0)-Ремг/(0) = 0, в'(0) - Рейв (0) = О, (6.3.40а)

/(!) = О, в'(1)-=0. (6.3.40Ь)

Уравнения (6.3.33) для этого случая (в компонентах v, w) записываются так:

vl-v+Ha

Рем 1

дЕ , дЕ

\.ду дЕ

pi<-t.;-fBDa[..+

д&

дЕ dQ

= -SVu (6.3.41) Wi =SOr,

- РйУд = -Si,

- Wi = SW.

Имеются в виду дополнительные уравнения типа (6.3.37). - Прим.

(6.3.42а)

Граничные условия для функций v п w получаются нз условий (6.3.40) в виде

v (0) - Pev (0) = О, v{ (0) - Pev, (0) = О,

w (0) - Реи^ (0) = О, w[ (0) - Pew, (0) = О,

v{l) = 0, v[{l) = 0,

Ч(1) = 0, w[il) = 0. ( -3-42b)

Для решения системы шести дифференциальных уравнений второго порядка (6.3.39), (6.3.41) воспользуемся методом стрель-бы. Недостающие начальные условия выберем в точке z = 1 (это связано с тем, что соответствующие задачи Коши неустой-чивы, если двигаться от точки 2 = 0 к точке 2=1, см. п. 6.1.2):

1/(1) = 11 6(1) = 112, ик(1) = Лз, Oi(l) = n4.

и^к(1) = т15, и^,(1) = Лб. (б--)

Вместе с условиями (6.3.40Ь) и (6.3.42Ь) мы получаем 12 начальных условий, необходимых для интегрирования уравнений (6.3.39) н (6.3.41). Поскольку собственная функция

u = UR-f iui = (0r, wji) + i{vi, Wi)

определена с точностью до ненулевого комплексного множителя, для функций Ur и Ui можно зафиксироввть два каких-либо конкретных значения; например, можно задать Or(1) = Oi(l) = 1 (т. е. положить в (6.3.43) 113 = 114= 1). В точке 2 = 0 мы получаем шесть нелинейных уравнений (6.3.40а) и (6.3.42а) относительно шести неизвестных цх, г[2, Ць, Цв, s и Da. Эту систему можно решать методом Ньютона, причем матрица Якоби вычисляется с помощью уравнений в вариациях, либо исходя из соответствующих разностных формул. Пример сходимости метода Ньютона к точке комплексной бифуркации представлен в табл. 6.11. При этом соответствующие итерации сходятся к точке, которая лежит на ветви решений с высокой выходной степенью конверсии у(1) (ср. качественно близкую диаграмму решений на рис. 6.8 для аналогичных значений параметров). В указанной точке стационарное решение (при убывании Da) теряет устойчивость, и от него отходит ветвь устойчивых периодических решений (см. п. 6.5.1). Описанный здесь прямой итерационный метод был рассмотрен на примере двух уравнений (4.3.7) (или, соответственно, (Р14-7), (Р14-8)). Метод легко обобщается на задачи с большим числом уравнений, но размерность возникающих при этом (конечномерных) задач, очевидно, возрастает.

Таблица 6.11. Пример сходимости метода Ньютона к точке комплексной бнфуркацнн для задачи 14 (у->оо, S = 12, р = 2, 0 = О, Рен = 2, Рем = 2, Le = 1).

6.4. МЕТОДЫ ДИНАМИЧЕСКОГО МОДЕЛИРОВАНИЯ ПАРАБОЛИЧЕСКИХ УРАВНЕНИЙ')

Для решения нелинейных дифференциальных уравнений с частными производными параболического типа существует целый ряд численных методов. В большинстве случаев они представляют собой различные варианты конечно-разностных методов, метода прямых и метода конечных элементов. Наряду со специальными работами имеется много учебников и монографий, посвященных численным методам решения дифференциальных уравнений с частными производными [6.23-6.28]. Ниже мы лишь кратко опишем основные конечно-разностные подходы и обсудим проблему их эффективной алгоритмизации. Все рассмотрение будет проведено на примере системы типа реакция-диффузия (4.3.7), т. е. для случая одной пространственной переменной.

В полуполосе 0г^1, 0<оо выберем равномерную сетку узловых точек {zu tj);

Zt = ih, i = 0, 1, n, h=l/n, ti = h, j = 0, 1, 2, T>0.

(6.4.1)

Имеется в виду численное моделирование дннамнки, т. е. численное решение задачи Кошн (точнее, задачи с начальными н краевыми условиями).-Ярим, ред.

+ -W iyi-г - Ч + yU) + Sl - (6-4.3b)

Формулы (6.4.3) включают в себя явную схему (ш = 0), чисто неявную схему {w = \) (обе с погрешностью аппроксимации порядка 0(т4-Л2)), а также классическую схему Кранка - Николсона с погрешностью аппроксимации порядка 0(т2--+ Л2) (йу = 1/2)>. В случае граничных условий 1 рода естественно положить при всех /

= х/+ = х; у1=у, у'п^=у. (6.4.4)

При этом формулы (6.4.3) записываются для г=1,2, п- 1. В случае граничных условий второго рода (ГУ2) уравнения (6.4.3) рассматриваются для г = 0, 1, п, а граничные условия (4.3.12) заменяются соотношениями

xL\ = x[\ 4+1= x/t; = ynX\=ylt\. (6.4.5)

Схема (6.4.3) является одношаговой по переменной /, а это значит, что новый слой / -f- 1 можно вычислить исходя из старого слоя у. Это значит, что значения /+ ° и g зависят только от xi, х'+\ у' и у^+К Явная схема {w = 0) представляет собой наиболее простой случай, поскольку значения неизвестных функций на новом слое (с верхним индексом ;+ 1) непосредственно выражаются через значения на старом слое (с верхним индексом /). При этом нелинейные члены и gl легко аппроксимируются подстановкой значений xi и yi в функции f я g. Вместе с тем, в случае явной схемы должно выполняться условие численной устойчивости, накладывающее определенные

Справедливость этих утверждений зависит от того, как определены fi+w и gi+r-Прим. ред.

Значения приближенного решения в этих узловых точках будем обозначать как

x{~x{z.,t.), y\~y{z t.). (6.4.2)

Наиболее часто используется шеститочечная разностная схема типа Кранка - Николсона. Для системы (4.3.7) она дает:

+ (1-1 - 22- + + /Г . (6.4.3а)

1 , . DvW

ограничения на величину шага т, а именно

< (w ж)- --

Эффективность указанной явной схемы может оказаться довольно низкой, если из условия (6.4.6) приходится выбирать очень малый шаг т.

При 1ЮфО соотношение (6.4.3) представляет собой уже систему уравнений относительно неизвестных с верхним индексом Такой метод называется неявным. При метод оказывается численно устойчивым, причем длина шага т ничем не ограничивается. Записав нелинейные члены /[+ , gl в уравнениях (6.4.3) в виде

fi. + -~f{wx{++{l-w)xl, wy>+ + {l-w)yl) (6.4.7)

(член записывается аналогично), мы получаем систему

2/г-f-2 нелинейных уравнений, которую можно решать итерациями (например, методом Ньютона). Поскольку у нас имеется хорошее начальное приближение (а именно значения на старом слое /), обыкновенно оказывается достаточным одной-двух итераций. Метод Ньютона предпочтительнее использовать при анализе быстрых химических реакций (т. е. тогда, когда производные f я g велики. - Ред.). При этом число итераций может оказаться более высоким. Для вычисления нового слоя используются также неитерационные подходы, в которых нелинейные члены аппроксимируются иначе - исходя из уже найденных значений неизвестных. По существу в таких подходах мы имеем дело с одной итерацией, проводимой каким-либо методом, ана-.логичным методу Ньютона.

Остановимся на этом вопросе несколько подробнее.

6.4.1. Замена нелинейных членов

Простейшая разностная замена нелинейных членов в схеме (6.4.3) заключается в использовании линеаризации по и т. е. для члена /!+° мы полагаем

fi+w ai + bixi+ + ciyi-\ (6.4.8)

где коэффициенты а(, Ы., с\ зависят от известных значений решения на старом (/-м) слое. Так, например, при ш=1 член