[5.8] Ralston А.: А First Course in Numerical Analysis, McGraw-Hill, New York, 1965.

[5.9

[5.io;

Kubicek M.: Appl. Math. Comput. 1 (1975). 341. Kubicek M., Marek M.: Computational Methods in Bifurcation Theory and Dissipative Structures, Springer, New York, 1983.

5.П] Kubicek M., Marek M.: Appi. Math. Comput. 5 (1979), 250. 5.12] Kubicek M., Klfc A.: Appi. Math. Comput. 13 (1983), 125. 5.13} Kubicek M., StuchI 1., Marek M.: J. Comput. Phys. 48 (1982), 106. 5.14] Kubicek M: Chem. Engng. Sci. 34 (1979), 1078. 5.15] Kubicek M.: SIAM J. Appl. Math. 38 (1980), 103. 5.16 Holodniok M., KubiSek M.: Appl. Math. Comput, 15 (1984), 261. ;5.17J Holodniok M., Kubicek M.: J. Comput. Phys., 55 (1984), 254. 5.18] Marek M., Schreiber I.: Chaotic Behaviour of Deterministic Dissipative

Systems. Academia, Praha, 1984, and Cambridge University Press,

Cambridge, 1991.

[5.19] Stoer J., Bulirsch R.: Introduction to Numerical Analysis. Springer, New York, 1980.

[5.20] Holodniok M., Kubicek M.: Continuation of Periodic Solutions in Ordinary Differential Equations with Application to the Hodgkin-Huxley Model. In: Ortiz E. L., Ed.: NumericaI Approximation of Partial Differential Equations*, Elsevier, 1987, 397.

5.2I] Henon M.: Physica D (1982), 412.

5.22] Shimada I., Nagashima Т.: Prog. Theor. Phys. 61 (1979) 1605. 5.23] Benettin G., Froeschle C, Scheidecker J. P.: Phys. Rev. A 19 (1979), 2454.

[5.24] Rheinboldt W. C, Burkardt J. V.: ACM Trans. Math. Software 15 (1983), 215.

[5.25] Keller H. В.: Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems. In: Rabinowitz P. H., Ed.: Applications of Bifurcation Theory, Academic Press, New York, 1977.

5.26 Kubicek M., Holodniok M.: J. of Comput. Physics 70 (1987) 203.

5.27] Feigenbaum M. J.: J. Statist. Physics 6 (1979), 669.

5.28] May R.: Nature 261 (1976), 459.

Preston C: Iterates of maps on an interval. Springer, Berlin, 1983.

[5.29] Schreiber I., Holodniok M., Kubicek M., Marek M.: J. of Stat. Physics. 43 (1986), 489.

[5.30] Holodniok M., Kubicek M.: Math, and Computers in Simulation, 29 (1987), 33.

5.3П Linniger W., Willoughby R. A.: SIAM J. Num. Anal. 6 (1969), 47.

5.32 Kubicek M., ViSnak K.: Chem. Engng. Commun. I (1974), 291.

5.33] Hassard B. D., Kazarinoff N. D., Wan Y. H.: Theory and Applications of Hopf Bifurcation Cambridge University Press, 1981. [Имеется перевод: Хэссард Б, Казарннов Н., Вэн И. Теория и приложения бифуркации рождения цикла. - М.: Мир, 1985. - 280 с] [5.34] Brenan К. Е., Campbell S. L., Petzold L. R.: Numerical Solution of Initial - Value Prob ems in Differential - Algebraic Equations. North - Holland, Elsevier, New York, 1989. [5.35] Ахиезер H. И., Глазман И. М.: Теория линейных операторов в Гильбертовом пространстве. - М.: Наука, 1966.

ДОПОЛНИТЕЛЬНАЯ ЛИТЕРАТУРА

Джакалья Г. Методы теории возмущений для нелинейных систем. - М.:

Наука, 1979. -319 с. Джозеф Д. Устойчивость движений жидкости.-М.: Мир, 1981. -638 с.

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

ЧИСЛЕННЫЕ МЕТОДЫ

И АЛГОРИТМЫ,

ИСПОЛЬЗУЕМЫЕ

ДЛЯ АНАЛИЗА СИСТЕМ

С РАСПРЕДЕЛЕННЫМИ

ПАРАМЕТРАМИ

х(0) = х, у{0) = д, (6.1.3а)

х{1) = х, у{\) = у. (6.1.3b)

x(0) = 0, г/(0) = 0, (6.1.4а)

х'(1) = 0, /(1) = 0. (6.1.4Ь)

OxoJ (0) + Рхол: (0) = Ухо, Оуо/ (0) + М (0) = Yyo, (6.1.5а)

axiJ(l) + PxiA;(l) = Yxb ауУ(1) + М(1) = у1. (6.1.5Ь)

Очень часто системы с распределенными параметрами преобразуют в системы с сосредоточенными параметрами. При этом обычно используется метод прямых (см. п. 6.4.5) совместно с методом Галеркина, методом коллокаций или каким-либо разностным методом высокой точности [6.1, 6.2]. Иногда для указанного преобразования используются и так называемые спектральные методы [6.3]. При уменьшении погрешности аппроксимации (например, при выборе более мелкого шага дискретизации) возрастает размерность получающихся систем с сосредоточенными параметрами. Тем не менее для некоторых типов задач даже весьма грубая аппроксимация дает удивительно хорошие (качественно правильные) результаты [6.4].

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

6.1. СТАЦИОНАРНЫЕ РЕШЕНИЯ (МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ)

Для иллюстрации методов нахождения стационарных решений параболических дифференциальных уравнений с частными производными рассмотрим системы типа реакция-диффузия (см. задачи И-13). Левые части уравнений (4.3.7) при этом мы полагаем равными нулю. Таким образом, стационарное решение удовлетворяет уравнениям ( = d/dz)

x + fix, У) = 0, (6.1.1)

y + gix. У) = 0. (6.1.2)

Запишем граничные условия (4.3.8), (4.3.12), (4.3.13) (для граничных условий первого рода мы рассмотрим только симметричный случай (4.3.9)):

ГУ1 ГУ2 ГУЗ

yi-i-2yi-\-yi+i-\--g(Xi, Уд = 0, i=l, 2.....n-1.

(6.1.7) - 1. (6.1.8)

Если ДЛЯ граничных условий первого рода f{x,y) = g{x,y) = (i, то имеется однородное по пространству решение системы (6.1.1-2): x{z)x, y{z)y. Это решение, очевидно, удовлетворяет и граничным условиям второго рода. Численными решениями нелинейных краевых задач занимался целый ряд авторов, среди публикаций которых можно найти самые разные работы - от чисто теоретических статей до сугубо прикладных исследований. Здесь мы рассмотрим указанную проблему сравнительно кратко; читателей же, которые заинтересуются этой проблемой более глубоко, мы отсылаем к монографической литературе [6.5, 6.6, 6.7, 6.8, 6.9, 6.10].

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

6.1.1. Разностные методы

На промежутке ze[0, 1] выберем сетку узловых точек 2о = 0, z\, Zn=l. Сетка, как правило, выбирается эквидистантной, Zi = ih, i = 0, 1, п, где h=\/n есть шаг сетки. На этой сетке узлов величины х,-, у,- мы рассматриваем как аппроксимации значений решения x{Zi), y{zi), г = 0, 1, п. Далее, производные в уравнениях (6.1.1) и (6.1.2) заменяются соответствующими разностными формулами, например, трехточечными центральными разностями:

о , (6.1.6)

Подстановка этих выражений в уравнения (6.1.1) и (6.1.2) приводит к следующей системе нелинейных уравнений:

Xt i-2Xi-\-Xt+i-\--f(Xi, уд = 0, i=l, 2.....п-1,

Добавляя к ним аппроксимации граничных условий, мы получаем систему 2п-\-2 уравнений. Так, для случая ГУ1 имеем

0 = . Уй = У, (6.1.9а)

Хп = х, Уг, = д. (6.1.9Ь)

Для ГУ2 (или ГУЗ) мы должны, кроме того, заменить производные в крайних точках. Если для аппроксимации первой производной в ГУ2 использовать двухточечную замену типа

то порядок аппроксимации в формулах (6.1.7), (6.1.8) понизится с 0{h) до 0(h). Одна из возможностей сохранить порядок аппроксимации равным 0(ft2) - использовать несимметричную разностную замену первой производной, включающую три узловые точки:

, -Зд;о f Axi - Хг

Другая возможность состоит в использовании аппроксимации, основанной на виртуальных точках с индексами i = - 1 и i = = п-\-1. Именно, распространим соотношения (6.1.7) и (6.1.8) на случай индексов i = 0 и i = n и заменим граничные условия (6.1.4) трехточечными центральными разностями:

o-4f = 0 о~Чг = 0 (6.1.10а)

Подставляя значения x i, г/ 1, Xn+i и уп+i из этих соотношений в формулы (6.1.7) и (6.1.8), находим

2xi-2xo + -f{xo, Уо) = 0, (6.1.11а)

2г/,-2г/о + (л;о, г/о) = 0, (6.1.lib)

2л; ,-2л; + /(л; , г/ ) = 0, (6.1.12а)

2уп-1-2уп-1-§{Хп, Уп) = 0. (6.1.12Ь)

Так же, как и в случае ГУ1, мы получили систему 2п-\-2 нелинейных уравнений (6.1.7), (6.1.8), (6.1.11), (6.1.12) относительно 2п-\-2 неизвестных. Упорядочим эти неизвестные

(6.1.14)

2п-1. (6.1.7) при i = n~\,

2п. (6.1.8) при г = п- 1,

2п+ 1. (6.1.12а),

2п + 2. (6.1.12Ь).

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

хххО.........О

ххОхО........О

хОхххО.......О

ОхххОхО......О

ООхОхххО.....О

ОООхххОхО. . . .0

Читатель может легко достроить матрицу размещения и убедиться, что она является пятидиагональной (это означает, что крестики располагаются на главной диагонали и на четырех соседних диагоналях). Отметим, что в данном случае матрица Якоби системы также будет пятидиагональной. Если теперь для

следующим образом:

X = (j:o. г/о. 1. Уь Х2, х„ г/ г/ ). (6.1.13)

Перепишем теперь соответствующие уравнения для случая ГУ2 в виде последовательности:

1. (6.1.11а),

2. (6.1.1 lb),

3. (6.1.7) при г= 1,

4. (6.1.8) при г = 1,

5. (6.1.7) при г = 2,

6. (6.1.8) при 1 = 2,

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

Применим теперь метод конечных разностей для нахождения стационарного решения задачи 16. Это стационарное решение описывается соотношениями (см. уравнения (Р16-11), (Р16-12))

У + f / - Ф'У ехр=0, (6.1.15а)

+ -f е- + у^Ф'у ехр , Д^ = О, (6.1..15Ь)

а граничные условия (Р16-13), (Р16-14) принимают вид

г/(0) = 0, (6.1.16а)

в'(0) = 0, (6.1.16Ь)

у{1)-\--у'И)=1, (6.1.16с)

®(1)+ЖГ®(1) = 0- (6.1.16d)

Иногда при вычислении стационарного решения удается понизить размерность исходной задачи. Продемонстрируем эту возможность на данном примере, предполагая Nu = Sh. Умножая уравнение (6.1.15а) на множитель у^, складывая результат с уравнением (6.1.15Ь) и вводя обозначение

м = ург/ + в, (6.1.17)

мы получаем уравнение

+ -7 = 0. (6.1.18)

Граничные условия (6.1.16а, Ь, с, d) переписываются в виде

ы'(0) = 0, (6.1.19а)

() + W () = yP- (6.1.19b)

Левая часть уравнения (6.1.18) равняется [г°и']/г°, откуда с использованием (6.1.19а) находим ы' = 0. Интегрируя, имеем ы = С, и из условия (6.1.19Ь) следует, что и{г) = vp. Теперь из формулы (6.1.17) вытекает зависимость между в и у

в = ур(1-г/). (6.1.20)

(1 + Р(1-г ))Ч^Р 1 + Р(1-г/,-)

yP(i-,)

i= 1, 2, n-1.

Заметим, что при Nu ф Sh это соотношение оказывается несправедливым. Таким образом, система (6.1.15) - (6.1.16) сводится к краевой задаче для одного дифференциального уравнения второго порядка

+ у'- Ф'у ехр TTFTT -2)

с граничными условиями (6.1.16а, с). Рассмотрим далее для простоты предельный случай Nucxs, Shoo (см. гл. 4, задача 16). Выберем сетку узловых точек следующим образом: п = ih, г = О, 1, ..., п, h = \/n. Тогда разностный аналог уравнения (6.1.21) принимает вид

ih Ih *i/xP 1 + р(1 ) и,

/=1, 2, п- 1. (6.1.22а)

Разностный аналог граничного условия (6.1.16а) получим, введя виртуальную точку г-\. В случае г = Ь, однако, в уравнении (6.1.21) имеется неопределенное выражение {а/г)у' типа О/О. При / -> О находим

у + {а1г)у'Ма-\)у .

Далее, используя разностную замену для уравнения (6.1.21) при / = О, с учетом требования у-i - у\ (которое вытекает из условия (6.1.16а)) получаем окончательно

{2у, - 2у,) - Ф^уоехр .ll] =0. (6.1.22b)

Наконец, граничное условие (6.1.16с) при Sh- -oo аппроксимируется с помощью соотношения

г/ -1=0. (6.1.22с)

Читатель может легко построить матрицу Якоби G = [§ ;/],

г, / = О, 1..... п, для решения системы методом Ньютона.

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

1 а . 10 1

gi,i-l --2 2ih - .....

i+i =-ХГ + -о^> i=U 2, П-\,

Соответствующие выражения для 00, 01, gnn. нетрудно получить дифференцированием уравнений (6.1.22Ь, с). Ход итерационного процесса в случае применения метода Ньютона представлен в табл. 6.1. Точное значение у(0) в этой задаче равно 0,5521, так что погрешность аппроксимации при /г = 0,1 сказывается лишь в четвертом знаке. Заметим, что приведенный выше пример оказывается настолько простым, что его вполне можно анализировать с помощью небольшой персональной ЭВМ.

В качестве более сложного примера рассмотрим решение задачи 17. Соответствующие разностные аналоги уравнений (Р17-16) -(Р17-20) имеют вид

Я„ = 0, Fn = 0, Gn=S.

т

X О

s о

о н S л

S о

О

н а> s

Яо = 0,

F, = 0, (6.1.23а) is Go=l,

Hi - Я, 1 + Л VR (Fi + Fi i)=0,

-±h/RH,{F,,-F, ,)--hRe{FI-GI-\-k) = 0, Gj i - 2G, + G,-+i - 2Л2 Re --h VR(G,+ 1 - G, i) = 0, 1 i=l, 2, n-l, (6.1.23b) 1 Hn - Я„ 1 + h VR(/ +/ ,)=0,

(6.1.23c)

s -

a.- с

8в о ю

(6.1.23d)

ее н

S о

,2 2.

Допустим, что значения параметров Re и S заданы. Упорядочим неизвестные следующим образом: Яо, F, Go, Hx, Fx, Gi, ..., Hn-x, Fn-i, Gn-x, Hn, Fn, Gn, k. Мы получим тогда семи-диагональную матрицу размещения с полностью заполненным последним столбцом (что соответствует появлению неизвестной k во всех уравнениях (6.1.23Ь)). В этом случае общее число уравнений оказывается равным 3(п-1) + 7. При интересных для этой задачи числах Рейнольдса (Re 500 1000) нужно довольно много узлов. Например, для Re = 625 необходимо л= 100-=-200 [6.11]. Подобную задачу следует решать уже на достаточно мощных ЭВМ, даже если использовать специальную программу для решения систем линейных алгебраических уравнений с почти семидиагональной матрицей. Кроме того, в данной задаче при достаточно больших значениях параметра Re (Re > 700) появляются паразитные решения. Так, при л = 100 и л = 200 существуют решения, которые в случае более мелкого разбиения, например при п - 800, исчезают (см. [6.11]). Учитывая это обстоятельство, при использовании разностных методов представляется необходимым результаты, полученные при достаточно грубом разбиении, пересчитывать на более мелкой сетке узловых точек. При этом только хорошее совпадение результатов для нескольких последовательных дроблений может служить критерием правильности решения.

6.1.2. Метод стрельбы

Основная идея метода стрельбы заключается в сведении краевой задачи к задаче Коши, для решения которой можно использовать стандартные программные средства, имеющиеся практически на каждой ЭВМ (см. § 5.7). Ниже мы рассмотрим простейшие варианты метода стрельбы, вполне достаточные для понимания его концептуальной стороны. Более подробное изложение этого метода читатель может найти в учебной литературе [6.7, 6.8].

Опишем подробнее метод стрельбы на примере задачи (6.1.1), (6.1.2) с граничными условиями (6.1.4). Для того чтобы данную задачу четвертого порядка можно было решать как задачу Коши, мы должны задать в некоторой точке четыре начальных условия. Воспользуемся тем, что в точке z = О уже имеются два заданных условия, а именно условия (6.1.4а), и выберем два дополнительных условия вида

т = Чь У(0) = %. (6.1.24)

Полученную задачу Коши на промежутке от z == О до z = 1 можно решать с помощью какого-либо из методов, рассмотрен-