Кирьянов Д. - MathCad 11 (1077323), страница 47
Текст из файла (страница 47)
Решение системы Лоренца с измененным параметром г=10Замечательно, что решение подобных нелинейных динамических системможно получить только численно, поэтому их изучение стало бурно развиваться с ростом возможностей вычислительной техники в последние полвека.11.4. Фазовый портретдинамической системыДо сих пор в этой главе в качестве примеров расчета динамических системмы приводили графики, траекторий на фазовой плоскости.
Однако для надежного исследования фазового портрета необходимо решить систему ОДУбольшое количество раз с самыми разными начальными условиямиЧасть III. Численные методы288(и, возможно, с разным набором параметров модели), чтобы посмотреть,к каким аттракторам сходятся различные траектории. В Mathcad можнореализовать эту задачу, например, в форме алгоритма, приведенного в листинге 11.10 для решения системы уравнений автокаталитической химической реакции с диффузией.
Эта модель, называемая моделью брюсселятора,предложена в 1968 г. Лефевром и Пригожиным. Неизвестные функции отражают динамику концентрации промежуточных продуктов некоторой реальной химической реакции. Параметр модели в равен исходной концентрации катализатора.Листинг 11.10. Построение фазового портрета для модели брюсселятрраоо2.51.50.5001v :=11.Е1 20.10.50.10.20.5 1.5В := 0 .
5D f t , у) :=t 0 := 0t l := 10М := 100U :=(о)У <г- VZ «- rkfixed (у, tO, tl,M,D)Z1(1)Z1Z1(2)f- Z(2)17 т \(1)for ke 1 .. last [\v )у ^ v (k)Z *- rkfixed (y,tO,tl,M,D)(0)(0)Z2 (1)(1)Z2Zls t a c k ( Zl , Z2ZlПредложенный алгоритм формирует из отдельных матриц решений системыОДУ с разными начальными условиями объединенную матрицу и. ПарыГлава 11.
Обыкновенные дифференциальные уравнения289начальных условий задаются в первой строке листинга в виде матрицы vразмера 2хю. Последнее означает построение десяти траекторий. Для тогочтобы поменять количество траекторий, измените соответствующим образом размер этой матрицы. Затем (рис. 11.14) элементы матрицы и выводятсяна график в виде отдельных точек.
Отсутствие соединения точек линиямиявляется недостатком алгоритма, но это плата за возможность представить вMathcad несложным образом сразу большое количество траекторий на фазовой плоскости.1.5ООнОООСОО0.52.5Рис. 11.14, Фазовый портрет брюсселятора при В=0 . 5 (листинг 11.10)Как видно из рис. 11.14, все траектории, вышедшие из разных точек, асимптотически стремятся к одному и тому же аттрактору ( 1 , 0 . 5 ) . Из теориидинамических систем нам известно, что такой аттрактор называется узлом(с узлом мы уже встречались в примерах разд. 11.1). Конечно, в общем случаепри анализе фазового портрета желательно "прощупать" большее число траекторий, задавая более широкий диапазон начальных условий.
Не исключено, что в других областях фазовой плоскости траектории будут сходитьсяк другим аттракторам.Эволюцию фазового портрета брюсселятора можно наблюдать, проводя расчеты с различным параметром в. При его увеличении узел будет сначалапостепенно смещаться в точку с координатами (1,в), пока не достигнетбифуркационного значения в=2. В этой точке происходит качественная перестройка портрета, выражающаяся в рождении предельного цикла.
Придальнейшем увеличении в происходит лишь количественное изменениеЧасть III. Численные методы290параметров этого цикла. Решение, полученное при в=2.5, показано нарис. 11.15.ПримечаниеЧтобы найти аттракторы динамической системы, как известно, нужно решитьсистему алгебраических уравнений, получающуюся из системы ОДУ заменойнулями их левых частей. Эти задачи также удобно решать средствами Mathcad(см. гл.
8). В частности, исследование зависимости фазового портрета от параметров системы ОДУ и поиск бифуркаций можно проводить методами продолжения (см. разд. "Метод продолжения по параметру" гл. 8).Рис. 11.15. Фазовый портрет брюсселятора при в=2 . 5Читатели, сталкивающиеся с расчетом динамических систем, несомненнооценят возможности Mathcad по построению фазовых портретов и исследованию бифуркаций. Возможно также, что они найдут лучшие программные решения этой задачи, чем алгоритм, предложенный в данном разделеавтором.11.5.
Жесткие системы ОДУДо сих пор мы имели дело с "хорошими" уравнениями, которые надежнорешались численными методами Рунге-Кутты. Однако имеется класс такназываемых •жестких (stiff) систем ОДУ, для которых стандартные методыпрактически неприменимы, поскольку их решение требует исключительномалого значения шага численного метода.
Некоторые из специальных алгоритмов, разработанных для этих систем, реализованы в Mathcad.Глава 11 .Обыкновенные дифференциальные уравнения29111.5.1. Что такое жесткие ОДУ?Строгого общепринятого математического определения жестких ОДУ нет.В рамках этой книги будем считать, что жесткие системы — это те уравнения, решение которых получить намного проще с помощью определенныхнеявных методов, чем с помощью явных методов (типа тех, что были рассмотрены в предыдущих разделах). Примерно такое определение былопредложено в 1950-х годах классиками в этой области Кертиссом и Хиршфельдером. Начнем разговор о жестких ОДУ с примера нежесткого уравнения (листинг 11.11), решение которого показано на рис.
11.16. На том жеграфике показано решение подобного ОДУ с другим коэффициентом приправой части, равным не -ю, а -зо. Решение обоих уравнений не составилотруда, и численный метод Рунге-Кутты дал правильный результат.| Листинг 11.11. Решение нежесткого ОДУ;] |Gi v e n—у (t)= - Ю - (у ( t ) - c o s ( t ) )dtУ(0)= 1у :- O d e s o l v e ( t , 1)— у ( t ) = - 1 0 .
(у ( t ) - c o s ( t ) )_ у ( t ) = - 3 0 - ( y ( t ) - c o s ( t ) Г-Рис. 1 1 . 1 6 . Решение нежестких ОДУ методом Рунге-Кутты(листинг 11.11)На рис. 11.17 показано решение того же ОДУ с коэффициентом -50. Вас,несомненно, должен насторожить результат, выданный Mathcad.
Характерная "разболтка" решения говорит о неустойчивости алгоритма. Первое, чтоЧасть III. Численные методы292можно сделать, — увеличить количество шагов в методе Рунге-Кутты. Для этого достаточно добавить третий параметр step в функциюodesolve(t,i<step). После нескольких экспериментов можно подобратьтакое значение step, которое будет обеспечивать устойчивость решения.Читатель может самостоятельно убедиться, что при step>20 "разболтка"пропадает, и решение становится похожим на графики, показанные нарис. 11.16.—у ( t )dtV (t)0.20.4=-5D. (у (t) - c o s (t))0.6Рис. 11.17. Неверное решение более жесткого ОДУ методом Рунге-КуттыТаким образом, во-первых, мы выяснили, что одни и те же уравнения сразными параметрами могут быть как жесткими, так и нежесткими. Вовторых, чем жестче уравнение, тем больше шагов в обычных численных методах требуется для его устойчивого решения.
С классическим примеромОДУ из листинга 11.11 все получилось хорошо, т. к. оно было не очень жестким, и небольшое увеличение числа шагов разрешило все проблемы. Длярешения обычными методами более жестких уравнений требуются миллионы, миллиарды и даже большее число шагов.ПримечаниеНекоторые ученые замечают, что в последние годы методы Рунге-Кутты сталиуступать свое главенствующее положение среди алгоритмов решения ОДУ методам, способным решать жесткие задачи.Исторически, интерес к жестким системам возник в середине XX века приизучении уравнений химической кинетики с одновременным присутствиемочень медленно и очень быстро протекающих химических реакций. Тогданеожиданно оказалось, что считавшиеся исключительно надежными методыГлава 11. Обыкновенные дифференциальные уравнения293Рунге-Кутты стали давать сбой при расчете этих задач.
Рассмотрим классическую модель взаимодействия трех веществ (Робертсон, 1966), которая какнельзя лучше передает смысл понятия жесткости ОДУ.Пусть вещество "О" медленно превращается в " i N : "0"-»"i" (со скоростьюo . i ) , вещество " 1 " при каталитическом воздействии самого себя превращается очень быстро в вещество "2": " i " + n i " - » " 2 " + " i " ( ю 3 ) . И, наконец, подобным образом (но со средней скоростью) реагируют вещества "2" и " 1 " :' Ч " + " 2 " ^ " 0 " + "2" ( ю 2 ) . Система ОДУ, описывающая динамику концентрации реагентов, с попыткой решения методом Рунге-Кутты, приведена влистинге 11.12.Листинг 11.12. Жесткая система ОДУ химической кинетики2-0.1.уо+ Ю -у1.у2F ( t , y ) :=D := rkf i x e d (yO , 0 , 50 , 20000 , F)Бросается в глаза сильно различающийся порядок коэффициентов при разных слагаемых.
Именно степень этого различия чаще всего и определяетжесткость системы ОДУ. В качестве соответствующей характеристики выбирают матрицу Якоби (якобиан) векторной функции F ( t , y ) , т . е . функциональную матрицу, составленную из производных F ( t , y ) (см. разд. "Частныепроизводные" гл. 7). Чем вырожденнее матрица Якоби, тем жестче системауравнений. В приведенном примере определитель якобиана и вовсе равеннулю при любых значениях у 0 , yi и у 2 (листинг 11.13, вторая строка).
В первой строке листинга 11.13 приведено напоминание способа вычисленияякобиана средствами Mathcad на примере определения элементов его первой строки.I Листинг 11.13. Якобиан рассматриваемой системы ОДУ химической кинетикиУо.1Эх-O.I10 •у20.1~102.у2-1030,310100 -у2—FЭх10 -уг- Ю02^-> о—FЭх-» 100-jЧасть ///. Численные методы294Для примера, приведенного в листинге 11.12, стандартным методом РунгеКутты все-таки удается найти решение (оно показано на рис. 11.18). Однакодля этого требуется очень большое число шагов, м=2оооо, что делает расчеты очень медленными.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
