Кирьянов Д. - MathCad 11 (1077323), страница 46
Текст из файла (страница 46)
Видно, что несмотря на высокую точность (10-5)и верный результат на конце интервала, левый график мало напоминаетправильный фазовый портрет (см. рис. 11.5 или правый график на рис. 11.6),начиная быть приемлемым только при предельно допустимом дляобсуждаемых функций значении асс=10-16.В заключение остановимся на влиянии выбора параметра асе на расчеты.Для этого воспользуемся простой программой, представленной на листинЮЗак.
984Часть III, Численные методы280ге 11.6. В ней из матрицы решения все той же задачи Коши взято лишь полученное значение одной из функций на правой фанице интервала. Но затоэтот результат оформлен в виде функции пользователя у (г), в качестве аргумента которой выбран параметр асе функции buistoer.0.1--0.1-0.050.1Рис. 11.6. Фазовый портрет, полученный b u i s t o e r (слева)и r k a d a p t (справа) (листинг 11.5)Листинг 11.6. Использование решения ОДУ для определенияфункции пользователяD ( t , у) :=у(е):=Y1уО :=-У0 - 0 . 10.1u«- buistoer (yO,0,50,E,D,30,O.OOOl)М <- l e n g t h (u '•*•'J- 1им. 1Вычисленный вид у(е) показан на рис.
11.7 вместе с аналогичным результатом для функции rkadapt. Как видно, в данном примере численные методыработают несколько по-разному. Метод Рунге-Кутты дает результат темближе к истинному, чем меньше выбирается е=асс. Метод Булирша-Штерадемонстрирует менее естественную зависимость у(е): даже при относительно больших е реальная точность остается хорошей (намного лучше методаРунге-Кутты).
Поэтому для экономии времени расчетов (подчеркнем ещераз: для данной конкретной задачи) в функции buistoer можно выбирать ибольшие асе.Чтобы обеспечить заданную точность, алгоритмы, реализованные во встроенных функциях, могут изменять как количество шагов, разбивающих ин-Глава 11. Обыкновенные дифференциальные уравнения281тервал ( t o , t i ) , так и их расположение вдоль интерв&та. Чтобы выяснить,на сколько шагов разбивался интервал при расчетах у(е)на рис. 11.7 длякаждого £, следует вычислить размер получающейся матрицы. Для этогоможно, например, определить подобные функции:м ( е ) : = l e n g t h ( r k a d a p t (yO , 0 , 50 , е , D , 1 0 0 , O . O O O l )У0U)11.008)1rkadapt(yO, о , 50 , £ , D , 3 0 , 0 .
0 0 0 l ) _ (0.0078 _-0.0076 -bulstoerlyO12 -ю" 400SO , eЭ, 3 0 , O.OOOl)-14•10-4б£"*Ёio~*o.ooiРис. 11.7. Зависимость расчетного значения одного из уравнений системы ОДУна конце интервала от параметра асе (листинг 11.6)j60-'§i11rkadapt(vO. 0 . 50"•"•€>-.,.D1000 00 Ol)- . . -..40rkadapt(vO. 0&4£102 1050 . e.D3 0 . о . 00 0 l )••*«—sulstoer(yO , 0 , 50 , e , D 3 0•1 ID610-ч8 •ID0 . 000l)0.
001ВРис. 11.8. Зависимость числа шагов от параметра асе численных методов282Часть ill. Численные методыСравнив два результата применения rkadapt для к=зо и ы о о , обратитевнимание (рис. 11.8), как еще один параметр — максимальное число шагов к,влияет на вид м(е). Заметим, что такие же изменения параметра к на расчетм(е) посредством функции buistoer влияют слабо.Таким образом, проводя тестовые расчеты для различных задач и подбираянаилучший набор параметров, можно существенно сэкономить ресурсыкомпьютера.
Конечно, проводить подобный анализ стоит в случаях, когдавремя расчетов играет важную роль.11.3.3. Некоторые примерыВ предыдущих разделах были использованы примеры исключительно линейных уравнений, т. е. содержащих только первую степень неизвестных функций и их производных. Между тем, многие нелинейные уравнения демонстрируют совершенно удивительные свойства, причем решение большинстваиз них можно получить лишь численно. Рассмотрим несколько наиболееизвестных классических примеров систем ОДУ, имея в виду, что читателюони могут пригодиться как в познавательных, так и в практических целях.Это модели динамики популяций (Вольтерры), генератора автоколебаний(Ван дер Поля), турбулентной конвекции (Лоренца) и химической реакциис диффузией (Пригожина).
Все они (впрочем, как и уже приведенные вышев этой главе) содержат производные по времени t и описывают динамикуразличных физических параметров. Задачи Коши для таких моделей называют динамическими системами, и для их изучения центральным моментомявляется анализ фазовых портретов, т. е. решений, получающихся при выборе всевозможных начальных условий.(ПримечаниеJВ большинстве примеров, изложенных ниже, для построения фазового портрета рассчитывается несколько решений для разных начальных условий.Ограничимся в дальнейшем минимальными комментариями и приведемлистинги и графики решений без подробного обсуждения.Модель "хищник—жертва"Модель взаимодействия "хищник—жертва" независимо предложили в 1925—1927 гг. Лотка и Вольтерра. Два дифференциальных уравнения (листинг 11.7)моделируют временную динамику численности двух биологических популяций жертвы у0 и хищника у1. Предполагается, что жертвы размножаются спостоянной скоростью с, а их численность убывает вследствие поеданияхищниками.
Хищники же размножаются со скоростью, пропорциональнойколичеству пиши (с коэффициентом г), и умирают естественным образомГлава 11. Обыкновенные дифференциальные уравнения283(смертность определяется константой D). В листинге рассчитываются трирешения D, G, P ДЛЯ разных начальных условий.Листинг 11.7. Модель "хищник—жертва"С := 0 . 1F ( t , y ) :=D := 1(r := 0.1С-у о -г-у 0 -у]_у - D - y ! + г-у о -У1М := 500tO := 0t l := 10010D := rkf ixed | [G ;- rkf ixed|, tO , tl , М , F8104P :- rkf ixed101.5, tO , tl , M , F, tO , tl , M , FМодель замечательна тем, что в такой системе наблюдаются циклическоеувеличение и уменьшение численности и хищника (рис. 11.9), и жертвы,так часто наблюдаемое в природе.
Фазовый портрет системы представляетсобой концентрические замкнутые кривые, окружающие одну стационарную точку, называемую центром. Как видно, модельные колебания численности обеих популяций существенно зависят от начальных условий — послекаждого периода колебаний система возвращается в ту же точку. Динамические системы с таким поведением называют негрубыми.Рис.
1 1 . 9 . График решения (слева) и фазовый портрет (справа)системы "хищник—жертва" (листинг 11.7)Часть III. Численные методы284АвтоколебанияРассмотрим решение уравнения Ван дер Поля, описывающего электрические колебания в замкнутом контуре, состоящем из соединенных последовательно конденсатора, индуктивности, нелинейного сопротивления и элементов,обеспечивающихподкачкуэнергииизвне(листинг 11.8).Неизвестная функция времени y{t> имеет смысл электрического тока, а впараметре ц заложены количественные соотношения между составляющимиэлектрической цепи, в том числе и нелинейной компонентой сопротивления.Листинг 11.8. Модель Ван дер Поля (ц-l)u ;= 1Given,2d:У f t !- ц - ( 1 - уdt(t)2).-dt+ У (t)= 0у{0) = 0 .
0 1у ' (0) = 0у : = 0 d e s o l v e (t , 30Рис. 11.10. График решения (слева) и фазовый портрет (справа)уравнения Ван дер Поля (листинг 11.8)Глава 11. Обыкновенные дифференциальные уравнения285Решением уравнения Ван дер Поля являются колебания, вид которых дляц.= 1 показан на рис.
11.10. Они называются автоколебаниями и принципиально отличаются от рассмотренных нами ранее (например, колебаний маятника в разд. 11.3.2) тем, что их характеристики (амплитуда, частота, спектр)не зависят от начальных условий, а определяются исключительно свойствами самой динамической системы. Через некоторое время расчетов послевыхода из начальной точки решение выходит на один и тот же цикл колебаний, называемый предельным циклом. Аттрактор типа предельного циклаявляется замкнутой кривой на фазовой плоскости. К нему асимптотическипритягиваются все окрестные траектории, выходящие из различных начальных точек, как изнутри (рис. 11.10), так и снаружи (рис.
11.11) предельногоцикла.Рис. 11.11. Решение уравнения Ван дер Поляпри других начальных условиях у=-2 , у ' =-3(ПримечаниеЕсли компьютер у Вас не самый мощный, то расчет фазового портрета срис. 11.10—11.11 в Mathcad может занять относительно продолжительное время, что связано с численным определением сначала решения y ( t ) , а потомего производной. Время расчетов можно было бы существенно сократить, еслииспользовать вместо вычислительного блока Given/Odesolve одну из встроенных функций, которые выдают решение в виде матрицы, например r k f i x e d .Аттрактор ЛоренцаОдна из самых знаменитых динамических систем предложена в 1963 г.
Лоренцем в качестве упрощенной модели конвективных турбулентных движе-Часть HI. Численные методы286ний жидкости в нагреваемом сосуде тороидальной формы. Система состоитиз трех ОДУ и имеет три параметра модели (листинг 11.9). Поскольку неизвестных функций три, то фазовый портрет системы должен определяться нена плоскости, а в трехмерном пространстве.I Листинг 11.9.
Модель Лоренцаа := 10г := 27b)31 10 ^уО :=10V1 0•7-У1-F t,y) :=• y 2 + J- - y o - y ito := 0t l := 3 0N = 1000D - rkfixed ( уО , t o , t lт := D, N , F)Y •=D_ Dn7 .=Z(X, V. 2)Рис. 1 1 . 1 2 . Аттрактор Лоренца (листинг 11.9)Решением системы Лоренца при определенном сочетании параметров(рис. 11.12) является странный аттрактор (или аттрактор Лоренца) —притягивающее множество траекторий на фазовом пространстве, которое повиду идентично случайному процессу. В некотором смысле, аттрактор Ло-Глава 11. Обыкновенные дифференциальные уравнения287ренца является стохастическими автоколебаниями, которые поддерживаются в динамической системе за счет внешнего источника.Решение в виде странного аттрактора появляется только при некоторых сочетаниях параметров. В качестве примера на рис. 11.13 приведен результатдля г=ю и тех же значений остальных параметров.
Как видно, аттракторомв этом случае является фокус. Перестройка типа фазового портрета происходит в области промежуточных г. Критическое сочетание параметров, прикоторых фазовый портрет системы качественно меняется, называется в теории динамических систем точкой бифуркации. Физический смысл бифуркации в модели Лоренца, согласно современным представлениям, описываетпереход ламинарного движения жидкости к турбулентному.х ;= D ^'T •- D <°>1z := D ^Y = D <2>14015 -X.20.---•—••»...10 -Yz(V " "02^С011102030т(X,Y,Z)Рис. 11.13.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
