Кирьянов. Самоучитель MathCad 11. Кирьянов д в

• Название: Кирьянов д в
• Анкор: Кирьянов. Самоучитель MathCad 11.pdf
• Дата: 28.04.2017
• Размер: 10.75 Mb.
• Имя файла: Кирьянов. Самоучитель MathCad 11.pdf
• Тип: Реферат; #6148
• Категория: Информатика. Вычислительная техника
• страница: 21 из 36

1 •
J
уО
0.1 0
-5
u buistoer (уО , 0 , 50 , 10 . D , 300 , 0 . 0001
:= length (u
- 1
= 21 50 7.638
\
= 7.638 2.648 х Ю- Чтобы попробовать альтернативный численный метод, достаточно в листинге ИМЯ фуНКЦИИ buistoer rkadapt.
Внимание!
Функции b u i s t o e r и rkadapt (те, что пишутся с маленькой буквы) не предназначены для нахождения решения в промежуточных точках интервала,
хотя они и выдают их в матрице-результате. На рис. 11.6 показаны фазовые портреты рассматриваемой системы ОДУ, полученные с помощью b u i s t o e результат листинга и с помощью rkadapt (при соответствующей замене третьей строки листинга Видно, что несмотря на высокую точность (и верный результат на конце интервала, левый график мало напоминает правильный фазовый портрет (см. рис. 11.5 или правый график на рис. начиная быть приемлемым только при предельно допустимом для обсуждаемых функций значении
В заключение остановимся на влиянии выбора параметра на расчеты.
Для этого воспользуемся простой программой, представленной на листин-

Часть III, Численные методы В ней из матрицы решения все той же задачи взято лишь полученное значение одной из функций на правой фанице интервала. Но зато этот результат оформлен в виде функции пользователя (г, в качестве аргумента которой выбран параметр функции -
0 . 1 11.6. Фазовый портрет, полученный b u i s t o e r (слева)
и r k a d a p t (справа) (листинг
Листинг
решения ОДУ для определения
функции
D , у) :=
- 0
уО :=
0 . 1
b u i s t o e r l e n g t h (u
- 1 Вычисленный вид у(е) показан на рис. 11.7 вместе с аналогичным результатом для функции rkadapt. Как видно, в данном примере численные методы работают несколько по-разному. Метод дает результат тем ближе к истинному, чем меньше выбирается Метод демонстрирует менее естественную зависимость у(е): даже при относительно больших реальная точность остается хорошей (намного лучше метода
Рунге-Кутты). Поэтому для экономии времени расчетов (подчеркнем еще раз для данной конкретной задачи) в функции buistoer можно выбирать и большие
Чтобы обеспечить заданную точность, алгоритмы, реализованные во встроенных функциях, могут изменять как количество шагов, разбивающих ин-

Глава Обыкновенные дифференциальные уравнения
тервал таки их расположение вдоль Чтобы выяснить,
на сколько шагов разбивался интервал при расчетах у(е)на рис. 11.7 для каждого следует вычислить размер получающейся матрицы. Для этого можно, например, определить подобные функции l e n g t h
, 0 , 50 , е , D , 1 0 0 , O . O O O l ) У 0 .
0 .
.008 0 2 4
0
•10 50 ,
£ ,
,
£
D , 3 0 ,
30
Ё
0 .
-
,
Зависимость расчетного значения одного из уравнений системы ОДУ
на конце интервала от параметра (листинг 0
-
0
j
2 10
,
0 .
0 0 ,
10 50 5 0 50
,
8
. D
. D
, D
100 3 0 0
. о 0 .
. 0 0 0 0 0 Рис Зависимость числа шагов от параметра численных методов

Часть ill. Численные методы
Сравнив два результата применения rkadapt для и обратите внимание (рис. 11.8), как еще один параметр — максимальное число шагов к,
влияет на вид Заметим, что такие же изменения параметра на расчет посредством функции влияют слабо.
Таким образом, проводя тестовые расчеты для различных задачи подбирая наилучший набор параметров, можно существенно сэкономить ресурсы компьютера. Конечно, проводить подобный анализ стоит в случаях, когда время расчетов играет важную роль Некоторые примеры
В предыдущих разделах были использованы примеры исключительно
ных уравнений, те. содержащих только первую степень неизвестных функций и их производных. Между тем, многие нелинейные уравнения демонстрируют совершенно удивительные причем решение большинства из них можно получить лишь численно. Рассмотрим несколько наиболее известных классических примеров систем ОДУ, имея ввиду, что читателю они могут пригодиться как в познавательных, таки практических целях.
Это модели динамики популяций (Вольтерры), генератора автоколебаний
(Ван дер Поля, турбулентной конвекции (Лоренца) и реакции с диффузией (Пригожина). Все они (впрочем, как и уже приведенные выше в этой главе) содержат производные повремени и описывают динамику различных физических параметров. Задачи Коши для таких моделей называют динамическими системами и для их изучения центральным моментом является анализ фазовых портретов, те. решений, получающихся при выборе всевозможных начальных условий Примечание
В большинстве примеров, изложенных ниже, для построения фазового портрета рассчитывается несколько решений для разных начальных условий.
Ограничимся в дальнейшем минимальными комментариями и приведем листинги и графики решений без подробного обсуждения.
Модель
Модель взаимодействия независимо предложили в гг. Лотка и Вольтерра. Два дифференциальных уравнения (листинг моделируют временную динамику численности двух биологических популяций жертвы и хищника Предполагается, что жертвы размножаются с постоянной скоростью с, а их численность убывает вследствие поедания хищниками. Хищники же размножаются со скоростью, пропорциональной количеству пиши (с коэффициентом и умирают естественным образом

Глава 11. Обыкновенные дифференциальные уравнения
283
(смертность определяется константой D). В листинге рассчитываются три решения D, G, P ДЛЯ разных начальных условий.
Листинг 11.7. Модель
С :=
D := 1
(
F ( t , y )
:= 500
tO := 0 t l := 100 10
D
rkf ixed | [ |,
, tl ,
, F
8
r := 0.1
G
rkf ixed
P :- rkf ixed
10 4
10 1.5
, tO , tl , , F
, tO , tl , M , Модель замечательна тем, что в такой системе наблюдаются циклическое увеличение и уменьшение численности и хищника (рис и жертвы,
так часто наблюдаемое в природе. Фазовый портрет системы представляет собой концентрические замкнутые кривые, окружающие одну стационарную точку, называемую центром Как видно, модельные колебания численности обеих популяций существенно зависят от начальных условий — после каждого периода колебаний система возвращается в туже точку. Динамические системы с таким поведением называют негрубыми.
Рис.
График решения (слева) и фазовый портрет (справа)
системы
(листинг 11.7)

Часть III. Численные методы
Автоколебания
Рассмотрим решение уравнения Ван дер Поля, описывающего электрические колебания в замкнутом контуре, состоящем из соединенных последовательно конденсатора, индуктивности, нелинейного сопротивления и элементов, обеспечивающих подкачку энергии извне (листинг
Неизвестная функция времени имеет смысл электрического тока, а в параметре заложены количественные соотношения между составляющими электрической цепи, в том числе и нелинейной компонентой сопротивле- ния.
Листинг
Модель Ван дер Поля 1
Given
,2
d t
= 0 . 0 1
(0) = у (t , 30
d t
+ У ( t ) = Рис. 11.10. График решения (слева) и фазовый портрет (справа)
уравнения Ван дер Поля (листинг

Глава Обыкновенные дифференциальные уравнения
285
Решением уравнения Ван дер Поля являются колебания, вид которых для 1 показан на рис. 11.10. Они называются автоколебаниями и принципиально отличаются от рассмотренных нами ранее (например, колебаний маятника в разд. 11.3.2) тем, что их характеристики (амплитуда, частота, спектр)
не зависят от начальных условий, а определяются исключительно свойствами самой динамической системы. Через некоторое время расчетов после выхода изначальной точки решение выходит на один и тот же цикл колебаний, называемый предельным циклом Аттрактор типа предельного цикла является замкнутой кривой на фазовой плоскости. К нему асимптотически притягиваются все окрестные траектории, выходящие из различных начальных точек, как изнутри (рис. 11.10), таки снаружи (рис предельного цикла Решение уравнения Ван дер Поля при других начальных условиях ,
( Примечание
Если компьютеру Вас не самый мощный, то расчет фазового портрета с рис в Mathcad может занять относительно продолжительное время, что связано счисленным определением сначала решения y ( t ) , а потом его производной. Время расчетов можно было бы существенно сократить, если использовать вместо вычислительного блока одну из встроенных функций, которые выдают решение в виде матрицы, например Аттрактор Лоренца
Одна из самых знаменитых динамических систем предложена в 1963 г. Ло- ренцем в качестве упрощенной модели конвективных турбулентных движе-

Часть HI. Численные методы
ний жидкости в нагреваемом сосуде тороидальной формы. Система состоит из трех ОДУ и имеет три параметра модели (листинг Поскольку неизвестных функций три, то фазовый портрет системы должен определяться не на плоскости, а в трехмерном пространстве Листинг 11.9.
Лоренца а 10
:= 27
уО
F
N
D
:=
10 10
t,y) :=
:= 0 1000
(
b
+
t l :=
,
,
3
-
3 0
t l , N , F т := D
Y
D
Z D
(X, V.
Аттрактор Лоренца (листинг
Решением системы Лоренца при определенном сочетании параметров
(рис. 11.12) является странный аттрактор (или аттрактор Лоренца) притягивающее множество траекторий на фазовом пространстве, которое по виду идентично случайному процессу. В некотором смысле, аттрактор Ло-

Глава Обыкновенные дифференциальные уравнения является стохастическими автоколебаниями, которые поддерживаются в динамической системе за счет внешнего источника.
Решение в виде странного аттрактора появляется только при некоторых сочетаниях параметров. В качестве примера на рис приведен результат для и тех же значений остальных параметров. Как видно, аттрактором в этом случае является фокус. Перестройка типа фазового портрета происходит в области промежуточных Критическое сочетание параметров, при которых фазовый портрет системы качественно меняется, называется в теории динамических систем точкой бифуркации Физический смысл бифуркации в модели Лоренца, согласно современным представлениям, описывает переход ламинарного движения жидкости к турбулентному 0 2 0 0
T •- D <°>
"
0 х ;= D
Y
D
15 -
10 -
20 т D
11.13. Решение системы Лоренца с измененным параметром
Замечательно, что решение подобных нелинейных динамических систем можно получить только численно, поэтому их изучение стало бурно развиваться с ростом возможностей вычислительной техники в последние полвека Фазовый портрет динамической системы
До сих пор в этой главе в качестве примеров расчета динамических систем мы приводили траекторий на фазовой плоскости. Однако для надежного исследования фазового портрета необходимо решить систему ОДУ
большое количество раз с самыми разными начальными условиями

Часть III. Численные методы
(и, возможно, с разным набором параметров модели, чтобы посмотреть,
к каким аттракторам сходятся различные траектории. В Mathcad можно реализовать эту задачу, например, в форме алгоритма, приведенного в листинге для решения системы уравнений автокаталитической химической реакции с диффузией. Эта модель, называемая моделью брюсселятора,
предложена в г. Лефевром и Пригожиным. Неизвестные функции отражают динамику концентрации промежуточных продуктов некоторой реальной химической реакции. Параметр модели в равен исходной концентрации катализатора.
Листинг 11.10. Построение фазового портрета для модели брюсселятрра
о о v
:=
0.5 В := 0 . 5
D , у) :=
2 . 5 1.5 0 0 0.5
t 0 := 0 100
:= 10 1 Е . 1 0.5 0 . 1 0.2
U :=
(о)
У
Z
rkfixed (у, tO,
Z1
Z1
Z1
Z
(2)
\
for
1 .. last
у rkfixed
Z2
Z2
Z l stack ( Zl , Предложенный алгоритм формирует из отдельных матриц решений системы
ОДУ с разными начальными условиями объединенную матрицу Пары

Глава 11. Обыкновенные дифференциальные уравнения
289
начальных условий задаются впервой строке листинга в виде матрицы размера Последнее означает построение десяти траекторий. Для того чтобы поменять количество траекторий, измените соответствующим образом размер этой матрицы. Затем (рис. 11.14) элементы матрицы выводятся на график в виде отдельных точек. Отсутствие соединения точек линиями является недостатком алгоритма, но это плата за возможность представить в несложным образом сразу большое количество траекторий на фазовой. ООО Фазовый портрет брюсселятора при В . 5 (листинг
Как видно из 11.14, все траектории, вышедшие из разных точек, асимптотически стремятся к одному и тому же аттрактору (1,0.5). Из теории динамических систем нам известно, что такой аттрактор называется узлом
(с узлом мы уже встречались в примерах разд Конечно, в общем случае при анализе фазового портрета желательно "прощупать" большее число траекторий, задавая более широкий диапазон начальных условий. Не исключено, что в других областях фазовой плоскости траектории будут сходиться к другим аттракторам.
Эволюцию фазового портрета брюсселятора можно наблюдать, проводя расчеты с различным параметром в. При его увеличении узел будет сначала постепенно смещаться в точку с координатами (в, пока не достигнет бифуркационного значения в. В этой точке происходит качественная перестройка портрета, выражающаяся в рождении предельного цикла. При дальнейшем увеличении в происходит лишь количественное изменение

Часть III. Численные методы
параметров этого цикла. Решение, полученное при показано на рис. 11.15.
Примечание
Чтобы найти аттракторы динамической системы, как известно, нужно решить систему алгебраических уравнений, получающуюся из системы ОДУ заменой нулями их левых частей. Эти задачи также удобно решать средствами см. гл. 8). В частности, исследование зависимости фазового портрета от параметров системы ОДУ и поиск бифуркаций можно проводить методами продолжения (см. разд. "Метод продолжения по параметру" гл. 8).
Фазовый портрет при . Читатели, сталкивающиеся с расчетом динамических систем, несомненно оценят возможности Mathcad по построению фазовых портретов и исследованию бифуркаций. Возможно также, что они найдут лучшие программные решения этой задачи, чем алгоритм, предложенный в данном разделе автором Жесткие системы ОДУ
До сих пор мы имели дело с "хорошими" уравнениями, которые надежно решались численными методами Однако имеется класс так называемых жестких (stiff) систем ОДУ, для которых стандартные методы практически неприменимы, поскольку их решение требует исключительно малого значения шага численного метода. Некоторые из специальных алгоритмов, разработанных для этих систем, реализованы в Mathcad.

Глава дифференциальные уравнения Что такое жесткие ОДУ?
Строгого общепринятого математического определения жестких ОДУ нет.
В рамках этой книги будем считать, что жесткие системы — это те уравнения, решение которых получить намного проще с помощью определенных неявных методов, чем с помощью явных методов (типа тех, что были рассмотрены в предыдущих разделах. Примерно такое определение было предложено в х годах классиками в этой области и Хирш- фельдером. Начнем разговор о жестких ОДУ с примера нежесткого уравнения (листинг 11.11), решение которого показано на рис На том же графике показано решение подобного ОДУ с другим коэффициентом приправой части, равным не а -зо. Решение обоих уравнений не составило труда, и численный метод Рунге-Кутты дал правильный результат Листинг Решение нежесткого ОДУ |
у ( t ) =
( у ( t )
( t ) )
d t
= у :-
( t ,
— у ( t )
(у ( t )
( t ) )
( y ( t ) - c o s ( t )
1 1 . 1 6 . Решение нежестких ОДУ методом Рунге-Кутты
На
11.17 показано решение того же ОДУ с коэффициентом -50. Вас,
несомненно, должен насторожить результат, выданный Характерная "разболтка" решения говорит о неустойчивости алгоритма. Первое, что

Часть Численные методы
можно сделать, — увеличить количество шагов в методе ты. Для этого достаточно добавить третий параметр step в функцию После нескольких экспериментов можно подобрать такое значение step, которое будет обеспечивать устойчивость решения.
Читатель может самостоятельно убедиться, что при step>20 "разболтка"
пропадает, и решение становится похожим на графики, показанные на у ( t )
( у ( t ) - c o s d t
0.2 0.4 0.6
Неверное решение более жесткого ОДУ методом
Таким образом, во-первых, мы выяснили, что одни и те же уравнения с разными параметрами могут быть как жесткими, таки нежесткими. Во- вторых, чем жестче уравнение, тем больше шагов в обычных численных методах требуется для его устойчивого решения. С классическим примером
ОДУ из листинга все получилось хорошо, т. к. оно было не очень жесткими небольшое увеличение числа шагов разрешило все проблемы. Для решения обычными методами более жестких уравнений требуются миллионы, миллиарды и даже большее число шагов.
Примечание
Некоторые ученые замечают, что в последние годы методы Рунге-Кутты стали уступать свое главенствующее положение среди алгоритмов решения ОДУ методам, способным решать жесткие задачи.
Исторически, интерес к жестким системам возник в середине XX века при изучении уравнений химической кинетики с одновременным присутствием очень медленно и очень быстропротекающих химических реакций. Тогда неожиданно оказалось, что считавшиеся исключительно надежными методы

Глава 11. Обыкновенные дифференциальные уравнения стали давать сбой при расчете этих задач. Рассмотрим классическую модель взаимодействия трех веществ (Робертсон, 1966), которая как нельзя лучше передает смысл понятия жесткости ОДУ.
Пусть вещество "О" медленно превращается в (со скоростью вещество " 1 " при каталитическом воздействии самого себя превращается очень быстро в вещество "2":
+
+
И, наконец,
добным образом (но со средней скоростью) реагируют вещества "2" и " 1 " :
Система ОДУ, описывающая динамику концентрации реагентов, с попыткой решения методом Рунге-Кутты, приведена в листинге Листинг Жесткая ОДУ химической кинетики :=
ixed (yO , 0 , 50 , 20000 , Бросается в глаза сильно различающийся порядок коэффициентов при разных слагаемых. Именно степень этого различия чаще всего и определяет жесткость системы ОДУ. В качестве соответствующей характеристики выбирают матрицу Якоби (якобиан) векторной функции F ( t , y ) , те. функциональную матрицу, составленную из производных (см. разд. "Частные
производные" гл. 7). Чем вырожденнее матрица Якоби, тем жестче система уравнений. В приведенном примере определитель якобиана и вовсе равен нулю при любых значениях и (листинг 11.13, вторая строка. Впервой строке листинга 11.13 приведено напоминание способа вычисления якобиана средствами на примере определения элементов его первой строки Листинг Якобиан рассматриваемой системы ОДУ химической кинетики Эх Эхо Часть ///. Численные методы
Для примера, приведенного в листинге стандартным методом Рунге- все-таки удается найти решение (оно показано на рис Однако для этого требуется очень большое число шагов что делает расчеты очень медленными. При меньшем числе шагов численному алгоритму не удается найти решение. В процессе работы алгоритма оно расходится, и вместо результата выдает ошибку о превышении предельно шого числа.
Еще один факт, на который стоит обратить внимание, — это различие в порядке величины получающегося решения. Как видно из рис. 11.18, концентрация первого реагента существенно (в тысячи раз) превышает концентрацию остальных. Это свойство также очень характерно для жестких систем.
Примечание
В принципе, можно было бы снизить жесткость системы "вручную, применяя масштабирование. Для этого нужно искусственно уменьшить искомую функцию к примеру, в тысячу раз, разделив все слагаемые в системе ОДУ, содержащие на 1000. После масштабирования для решения полученной системы методом Рунге-Кутты будет достаточно взять всего М шагов.