8.1. ОДУ первого порядка

Дифференциальное уравнение первого порядка может по определению содержать, помимо самой искомой функции у (t), только ее первую производную y'(t). В подавляющем большинстве случаев дифференциальное уравнение можно записать в стандартной форме (форме Коши):

у (t)=f(y(t),t), (1)

и только с такой формой умеет работать вычислительный процессор Mathcad. Правильная с математической точки зрения постановка соответствующей задачи Коши для ОДУ первого порядка должна, помимо самого уравнения, содержать одно начальное условие — значение функции y(t0) в некоторой точке t0. Требуется явно определить функцию y(t) на интервале от t0 до tx. По характеру постановки задачи Коши называют еще задачами с начальными условиями (initial value problem), в отличие от краевых задач.

Для численного интегрирования одного ОДУ у пользователя Mathcad 11 (начиная с версии Mathcad 2000 Pro) имеется выбор — либо использовать вычислительный блок Given/odesoive, либо встроенные функции, как в прежних версиях Mathcad. Первый путь предпочтительнее из соображений наглядности представления задачи и результатов, а второй дает пользователю больше рычагов воздействия на параметры численного метода. Рассмотрим последовательно оба варианта решения.


8.1.1. Вычислительный блок Given/Odesolve

Вычислительный блок для решения одного ОДУ, реализующий численный метод Рунге-Кутты, состоит из трех частей:

  • Given — ключевое слово;
  • ОДУ и начальное условие, записанное с помощью логических операторов, причем начальное условие должно быть в форме у (t1) = b;
  • odesoive(t, t1) — встроенная функция для решения ОДУ относительно переменной t на интервале (t0,t1).

Допустимо, и даже часто предпочтительнее, задание функции Odesolve (t, t1, step) с тремя параметрами, где step— внутренний параметр численного метода, определяющий количество шагов, в которых метод Рунге - Купы, будет рассчитывать решение дифференциального уравнения. Чем больше step, тем с лучшей точностью будет получен результат, но тем больше времени будет затрачено на его поиск. Помните, что подбором этого параметра можно заметно (в несколько раз) ускорить расчеты без существенного ухудшения их точности.

Пример решения задачи Коши для ОДУ первого порядка у' =у-у2 посредством вычислительного блока приведен в листинге 11.1.

Листинг 11.1. Решение задачи Коши для ОДУ первого порядка

Не забывайте о том, что вставлять логические операторы следует при помощи панели инструментов Boolean (Булевы операторы). При вводе с клавиатуры помните, что логическому знаку равенства соответствует сочетание клавиш <Ctrl>+<=>. Символ производной можно ввести как средствами панели Calculus (Вычисления), как это сделано в листинге 11.1, так и в виде штриха, набрав его с помощью сочетания клавиш <Ctrl>+<F7> (соответствующий пример будет приведен ниже в листинге 11.3.) Выбирайте тот или иной способ представления производной из соображений наглядности представления результатов — на ход расчетов он не влияет.

Mathcad требует, чтобы конечная точка интегрирования ОДУ лежала правее начальной: t0<t1 (в листинге 11.1 t0=0,t1=10), иначе будет выдано сообщение об ошибке. Как можно заметить, результатом применения блока Given/odesoive является функция y(t), определенная на промежутке (t0,t1). Следует воспользоваться обычными средствами Mathcad, чтобы построить ее график или получить значение функции в какой-либо точке указанного интервала, например: у(3)=0.691.

Пользователь имеет возможность выбирать между двумя модификациями численного метода Рунге-Кутты. Для смены метода необходимо нажатием правой кнопки мыши на области функции odesolve вызвать контекстное меню и выбрать в нем один из двух пунктов: Fixed (Фиксированный шаг) или Adaptive (Адаптивный). По умолчанию применяется первый из них, т. е. метод Рунге - Кутты с фиксированным шагом. Подробнее о различии этих методов сказано в разд. 11.3.


8.1.2. Встроенные функции rkfixed, Rkadapt, Bulstoer

Альтернативный метод решения ОДУ перешел из прежних версий Mathcad. Он заключается в использовании одной из встроенных функций rkfixed, Rkadapt или Bulstoer. Этот способ несколько проигрывает первому и в простоте, и в наглядности Поэтому я советую предпочесть вычислительный блок Given/odesoive Однако иногда приходится решать ОДУ первого порядка с помощью второго способа, например, при следующих обстоятельствах:

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

Поскольку решение вторым способом одного ОДУ мало чем отличается от решения систем ОДУ (см разд 11.3), приведем пример его использования в задаче из листинга 11.1. практически без комментариев (см. листинг 11.2) и с помощью одной из трех существующих для этих целей встроенных функций rkfixed. Обратите внимание только на необходимость явного задания количества точек интегрирования ОДУ м=100 в третьей строке листинга, а также на получение результата, в отличие от вычислительного блока, не в виде функции, а в виде матрицы размерности M X 2. Она состоит из двух столбцов в одном находятся значения аргумента t (от t0 до t1 включительно), а в другом соответствующие значения искомой функции y(t).

Листинг 11.2. Решение задачи Коши для ОДУ первого порядка вторым способом

В листинге 11.2. приведен пример не лучшего стиля Mathcad-программиро-вания Сначала переменной у присвоено значение скаляра у=0 1, а затем этой же переменной присвоено матричное значение (результат решения ОДУ) Старайтесь избегать такого стиля, который ухудшает читаемость программы и может приводить, в более сложных случаях, к трудно опознаваемым ошибкам Неплохим решением было бы назвать результат по-другому, например u.

График решения рассматриваемого уравнения показан на рис. 11.1. Обратите внимание, что он соответствует получению решения в матричном виде (листинг 11.2), поэтому по осям отложены соответствующие столбцы, выделенные из матрицы у оператором <>.

Пример, решенный в листингах 11.1—11.2, взят из области математической экологии и описывает динамику популяций с внутривидовой конкуренцией Сначала происходит рост численности популяции, близкий к экспоненциальному, а затем выход на стационарное состояние

Рис. 11.1. Решение уравнения y' =y-y2 (листинг 11.2)


8.2. ОДУ высшего порядка

Обыкновенное дифференциальное уравнение с неизвестной функцией y(t), в которое входят производные этой функции вплоть до y(N) (t), называется ОДУ N-ГО порядка. Если имеется такое уравнение, то для корректной постановки задачи Коши требуется задать N начальных условий на саму функцию y(t) и ее производные от первого до (N-1) го порядка включительно. В Mathcad 11 можно решать ОДУ высших порядков как с помощью вычислительного блока Given/odesolve, так и путем сведения их к системам уравнений первого порядка.

Внутри вычислительного блока:

  • ОДУ должно быть линейно относительно старшей производной, т. е. фактически должно быть поставлено в стандартной форме;
  • начальные условия должны иметь форму y(t)=b или y(N) (t)=b, а не более сложную (как, например, встречающаяся в некоторых математических приложениях форма у (t) +у' (t) = b)

В остальном, решение ОДУ высших порядков ничем не отличается от решения уравнений первого порядка (см. разд. 11.1), что иллюстрируется листингом 11.3. Как Вы помните, допустимо написание производной как в виде знака дифференциала (так в листинге 11.3 введено само уравнение), так и с помощью штриха (так введено начальное условие для первой производной). Не забывайте пользоваться булевыми операторами при вводе уравнений и начальных условий. Полученное решение y(t) показано на рис. 11.2.

Листинг 11.3. Решение задачи Коши для ОДУ второго порядка

Рис. 11.2. Решение уравнения осциллятора (листинг 11.3)

В листинге 11.3 решено уравнение затухающего гармонического осциллятора, которое описывает, например, колебания маятника. Для модели маятника y(t) описывает изменения угла его отклонения от вертикали, y'(t) — угловую скорость маятника, y"(t) — ускорение, а начальные условия, соответственно, начальное отклонение маятника у (0) =0.1 и начальную скорость у' (0)= 0.

Второй способ решения ОДУ высшего порядка связан со сведением его к эквивалентной системе ОДУ первого порядка. Покажем на том же примере из листинга 11.3, как это делается. Действительно, если формально обозначить y0(t)sy(t), a yi(t)sy'(t)=y0'(t), то исходное уравнение запишется через функции y0(t) и y1(t) в виде системы двух ОДУ:

Именно эта система решается в качестве примера в разд. 11.3. Таким образом, любое ОДУ N-ГО порядка, линейное относительно высшей производной, можно свести к эквивалентной системе N дифференциальных уравнений.

8.3. Системы ОДУ первого порядка

Mathcad требует, чтобы система дифференциальных уравнений была представлена в стандартной форме.

Задание системы эквивалентно следующему векторному представлению, где Y и У ' — соответствующие неизвестные векторные функции переменной t размера NXI, ар — векторная функция того же размера и (N+i) количества переменных (N компонент вектора и, возможно, t). Именно векторное представление используется для ввода системы ОДУ в среде Mathcad.

Для того чтобы определить задачу Коши для системы ОДУ, следует определить еще N начальных условий, задающих значение каждой из функций yi(t0) в начальной точке интегрирования системы t0. В векторной форме они могут быть записаны, где в — вектор начальных условий размера NXI, составленный из y1(t0).

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

Обратите внимание на необходимость векторной записи как самого уравнения, так и начального условия. В случае одного ОДУ соответствующие векторы имеют только один элемент, а в случае системы N>I уравнений — N.


8.3.1. Встроенные функции для решения систем ОДУ

В Mathcad 11 имеются три встроенные функции, которые позволяют решать поставленную в форме (2—3) задачу Коши различными численными методами.

  • rkfixed(y0, t0, t1, M, D) — метод Рунге-Кутты с фиксированным шагом,
  • Rkadapt(y0, t0, t1, M, D) — метод Рунге-Кутты с переменным шагом;
  • Buistoer(y0, t0, t1, M, D) — метод Булирша-Штера;
    • у0 — вектор начальных значений в точке to размера NXI;
    • t0 — начальная точка расчета,
    • t1 — конечная точка расчета,
    • M — число шагов, на которых численный метод находит решение;
    • D — векторная функция размера NXI двух аргументов — скалярного t и векторного у При этом у — искомая векторная функция аргумента t того же размера NXI.

Соблюдайте регистр первой буквы рассматриваемых функций, поскольку это влияет на выбор алгоритма счета, в отличие от многих других встроенных функций Mathcad, например Find=fmd (см разд 11.3.2)

Каждая из приведенных функций выдает решение в виде матрицы размера (M+1)х(N+1). В ее левом столбце находятся значения аргумента t, делящие интервал на равномерные шаги, а в остальных N столбцах — значения искомых функций y0(t) ,y1(t), .. ,yN-1(t), рассчитанные для этих значений аргумента Поскольку всего точек (помимо начальной) м, то строк в матрице решения будет всего M+1

В подавляющем большинстве случаев достаточно использовать первую функцию rkfixed, как это показано в листинге 11.4 на примере решения системы ОДУ модели осциллятора с затуханием (см. разд. 11 2).

Листинг 11.4. Решение системы двух ОДУ

Самая важная — это первая строка листинга, в которой, собственно, определяется система ОДУ. Сравните рассматриваемую систему (разд. 11.2.1), записанную в стандартной форме, с формальной ее записью в Mathcad, чтобы не делать впоследствии ошибок. Во-первых, функция D, входящая в число параметров встроенных функций для решения ОДУ, должна быть функцией обязательно двух аргументов. Во-вторых, второй ее аргумент должен быть вектором того же размера, что и сама функция D. В-третьих, точно такой же размер должен быть и у вектора начальных значений уо (он определен во второй строке листинга).

Не забывайте, что векторную функцию D(t,y) следует определять через компоненты вектора у с помощью кнопки нижнего индекса (Subscript) с наборной панели Calculator (Калькулятор) или нажатием клавиши <[>. В третьей строке листинга определено число шагов, на которых рассчитывается решение, а его последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (о, 50).

Как выглядит все решение, показано на рис. 11.3. Размер полученной матрицы будет равен (M+DX(N+I), т.е. Ю1хз. Просмотреть все компоненты матрицы и, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки. Как нетрудно сообразить, на этом рисунке отмечено выделением расчетное значение первого искомого вектора у0 на 12-м шаге u12.1=0.01. Это соответствует, с математической точки зрения, найденному значению У0(6.о)=0.07. Для вывода элементов решения в последней точке интервала используйте выражения типа Uм1=7.523x103.

Рис. 11.3. Матрица решений системы уравнений (листинг 11.4)

Обратите внимание на некоторое разночтение в обозначении индексов вектора начальных условий и матрицы решения В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором столбце — первой компоненты и т. д.

Чтобы построить график решения, надо отложить соответствующие компоненты матрицы решения по координатным осям" значения аргумента и<0> — вдоль оси х, а u<1> и u<2> — вдоль оси У (рис. 11.4). Как известно, решения обыкновенных дифференциальных уравнений часто удобнее изображать не в таком виде, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций. При этом аргумент входит в них лишь параметрически. В рассматриваемом случае двух ОДУ такой график — фазовый портрет системы — является кривой на фазовой плоскости и поэтому особенно нагляден. Он изображен на рис. 11.5 (слева), и можно заметить, что для его построения потребовалось лишь поменять метки осей на u<1> и u<2>, соответственно.

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

Рис. 11.4. График решения системы ОДУ (11.2—1) (листинг 11.4)

Рис. 11.5. Фазовый портрет решения системы ОДУ при M=100 (слева) и М=200 (справа) (листинг 11.4)

В общем случае, если система состоит из N ОДУ, то фазовое пространство является N-мерным. При N>3 наглядность теряется, и для визуализации фазового портрета приходится строить его различные проекции.

На том же рис. 11.5, справа, показан для сравнения результат расчета фазового портрета с большим числом шагов. Видно, что в этом случае обеспечивается лучшая точность и, в результате, решение получается более гладким. Конечно, при этом увеличивается и время расчетов.

При решении систем ОДУ многие проблемы могут быть устранены при помощи простой попытки увеличить число шагов численного метода В частности, сделайте так при неожиданном возникновении ошибки "Found a number with a magnitude greater than 10^307" (Найдено число, превышающее значение 10307) Данная ошибка может означать не то, что решение в действительности расходится, а просто недостаток шагов для корректной работы численного алгоритма

В заключение следует сказать несколько слов об особенностях различных численных методов. Все они основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации, получаются алгоритмы различной точности и быстродействия. В Mathcad использован наиболее популярный алгоритм Рунге-Кутты четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо rkfixed другие функции, благо сделать это очень просто, благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе.

Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо, либо существуют участки медленных и быстрых его изменений. Метод Рунге-Кутты с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений — частыми. В результате, для достижения одинаковой точности требуется меньшее число шагов, чем для rkf ixed. Метод Булирша-Штера Buistoer часто оказывается более эффективным для поиска гладких решений.


8.3.2. Решение систем ОДУ в одной заданной точке

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

Такая задача требует меньше ресурсов компьютера, чем решение системы ОДУ на всем интервале, поэтому в Mathcad имеются модификации встроенных функций rkadapt и Bulstoer. Они имеют несколько другой набор параметров и работают быстрее своих аналогов.

  • rkadapt(y0,t0,t1,асе, D, k, s) — метод Рунге-Кутты с переменным шагом;
  • buistoer(y0, t0, t1, acc, D, k, s) —метод Булирша-Штера;
    • у0 — вектор начальных значений в точке t0;
    • t0,t1 — начальная и конечная точки расчета;
    • асе — погрешность вычисления (чем она меньше, тем с лучшей точностью будет найдено решение; рекомендуется выбирать значения погрешности в районе 0.001);
    • D — векторная функция, задающая систему ОДУ;
    • k — максимальное число шагов, на которых численный метод будет находить решение;
    • s — минимально допустимая величина шага.

Как легко заметить, вместо числа шагов на интервале интегрирования ОДУ, в этих функциях необходимо задать точность расчета численным методом значения функций в последней точке. В этом смысле параметр асе похож на константу TOL, которая влияет на большинство встроенных численных алгоритмов Mathcad. Количество шагов и их расположение определяется численным методом автоматически, чтобы обеспечить эту точность. Два последних параметра нужны для того, чтобы пользователь мог искусственно повлиять на разбиение интервала на шаги. Параметр k служит для того, чтобы шагов не было чрезмерно много, причем, нельзя сделать k>1000. Параметр s — для того чтобы ни один шаг не был слишком малым для появления больших погрешностей при разностной аппроксимации дифференциальных уравнений внутри алгоритма. Эти параметры следует задавать явно, исходя из свойств конкретной системы ОДУ. Как правило, проведя ряд тестовых расчетов, можно подобрать их оптимальный набор для каждого конкретного случая.

Пример использования функции buistoer для того же примера (11.2—1) приведен в листинге 11.5. В его первых двух строках, как обычно, определяется система уравнений и начальные условия; в следующей строке матрице и присваивается решение, полученное с помощью buistoer. Структура этой матрицы в точности такая же, как и в случае решения системы ОДУ посредством уже рассмотренных нами ранее встроенных

функций (см. разд. 11.3.1). Однако в данном случае нам интересна только последняя точка интервала. Поскольку сделанное численным методом количество шагов, т. е. размер матрицы и, заранее неизвестно, то его необходимо предварительно определить. Это сделано в четвертой строке листинга, присваивающей это число переменной м, в этой же строке оно выведено на экран. В предпоследней строке листинга осуществлен вывод решения системы ОДУ на конце интервала, т. е. в точке t=s0 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующим местом вектора из предыдущей строки).

Листинг 11.5. Решение системы двух ОДУ

Чтобы попробовать альтернативный численный метод, достаточно в листинге 11.5 заменить имя функции buistoer на rkadapt.

Функции buistoer и rkadapt (те, что пишутся с маленькой буквы) не предназначены для нахождения решения в промежуточных точках интервала, хотя они и выдают их в матрице-результате. На рис. 11.6 показаны фазовые портреты рассматриваемой системы ОДУ, полученные с помощью buistoer (результат листинга 11.5) и с помощью rkadapt (при соответствующей замене третьей строки листинга 11.5). Видно, что несмотря на высокую точность (10-5) и верный результат на конце интервала, левый график мало напоминает правильный фазовый портрет (см. рис. 11.5 или правый график на рис. 11.6), начиная быть приемлемым только при предельно допустимом для обсуждаемых функций значении асс=10-16.

В заключение остановимся на влиянии выбора параметра асе на расчеты. Для этого воспользуемся простой программой, представленной на листинre 11.6. В ней из матрицы решения все той же задачи Коши взято лишь полученное значение одной из функций на правой границе интервала. Но зато этот результат оформлен в виде функции пользователя у(е), в качестве аргумента которой выбран параметр асе функции bulstoer.

Рис. 11.6. Фазовый портрет, полученный bulstoer (слева) и rkadapt (справа) (листинг 11.5)

Листинг 11.6. Использование решения ОДУ для определения функции пользователя

Вычисленный вид у(е) показан на рис. 11.7 вместе с аналогичным результатом для функции rkadapt. Как видно, в данном примере численные методы работают несколько по-разному. Метод Рунге-Кутты дает результат тем ближе к истинному, чем меньше выбирается е=асс. Метод Булирша-Штера демонстрирует менее естественную зависимость у (Б): даже при относительно больших е реальная точность остается хорошей (намного лучше метода Рунге-Кутты). Поэтому для экономии времени расчетов (подчеркнем еще раз: для данной конкретной задачи) в функции bulstoer можно выбирать и большие асе.

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

Рис. 11.7. Зависимость расчетного значения одного из уравнений системы ОДУ на конце интервала от параметра асе (листинг 11.6)

Рис. 11.8. Зависимость числа шагов от параметра асе численных методов

Сравнив два результата применения rkadapt для k=30 и k=100, обратите внимание (рис. 11.8), как еще один параметр — максимальное число шагов k, влияет на вид м(е). Заметим, что такие же изменения параметра k на расчет м(е) посредством функции bulstoer влияют слабо.

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


8.3.3. Некоторые примеры

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

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

Ограничимся в дальнейшем минимальными комментариями и приведем листинги и графики решений без подробного обсуждения.

Модель "хищник—жертва"

Модель взаимодействия "хищник—жертва" независимо предложили в 1925— 1927 гг. Лотка и Вольтерра. Два дифференциальных уравнения (листинг 11.7) моделируют временную динамику численности двух биологических популяций жертвы Y0 и хищника Y1. Предполагается, что жертвы размножаются с постоянной скоростью с, а их численность убывает вследствие поедания хищниками. Хищники же размножаются со скоростью, пропорциональной количеству пищи (с коэффициентом r), и умирают естественным образом (смертность определяется константой D). В листинге рассчитываются три решения D, G, р для разных начальных условий.

Листинг 11.7. Модель "хищник-жертва"

Модель замечательна тем, что в такой системе наблюдаются циклическое увеличение и уменьшение численности и хищника (рис. 11.9), и жертвы, так часто наблюдаемое в природе. Фазовый портрет системы представляет собой концентрические замкнутые кривые, окружающие одну стационарную точку, называемую центром. Как видно, модельные колебания численности обеих популяций существенно зависят от начальных условий — после каждого периода колебаний система возвращается в ту же точку. Динамические системы с таким поведением называют негрубыми.

Рис. 11.9. График решения (слева) и фазовый портрет (справа) системы "хищник—жертва" (листинг 11.7)

Автоколебания

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

Листинг 11.8. Модель Ван дер Поля (м=1)

Рис. 11.10. График решения (слева) и фазовый портрет (справа) уравнения Ван дер Поля (листинг 11.8)

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

Рис. 11.11. Решение уравнения Ван дер Поля при других начальных условиях у=-2, у =-3

Если компьютер у Вас не самый мощный, то расчет фазового портрета с рис. 11.10—11.11 в Mathcad может занять относительно продолжительное время, что связано с численным определением сначала решения y(t), а потом его производной. Время расчетов можно было бы существенно сократить, если использовать вместо вычислительного блока Given/Odesolve одну из встроенных функций, которые выдают решение в виде матрицы, например rkfixed.

Аттрактор Лоренца

Одна из самых знаменитых динамических систем предложена в 1963 г. Лоренцем в качестве упрощенной модели конвективных турбулентных движений жидкости в нагреваемом сосуде тороидальной формы. Система состоит из трех ОДУ и имеет три параметра модели (листинг 11.9). Поскольку неизвестных функций три, то фазовый портрет системы должен определяться не на плоскости, а в трехмерном пространстве.

Листинг 11.9. Модель Лоренца

Рис. 11.12. Аттрактор Лоренца (листинг 11.9)

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

Решение в виде странного аттрактора появляется только при некоторых сочетаниях параметров. В качестве примера на рис. 11.13 приведен результат для г=ю и тех же значений остальных параметров. Как видно, аттрактором в этом случае является фокус. Перестройка типа фазового портрета происходит в области промежуточных г. Критическое сочетание параметров, при которых фазовый портрет системы качественно меняется, называется в теории динамических систем точкой бифуркации. Физический смысл бифуркации в модели Лоренца, согласно современным представлениям, описывает переход ламинарного движения жидкости к турбулентному.

Рис. 11.13. Решение системы Лоренца с измененным параметром г=10

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


8.4. Фазовый портрет динамической системы

До сих пор в этой главе в качестве примеров расчета динамических систем мы приводили графики траекторий на фазовой плоскости. Однако для надежного исследования фазового портрета необходимо решить систему ОДУ большое количество раз с самыми разными начальными условиями (и, возможно, с разным набором параметров модели), чтобы посмотреть, к каким аттракторам сходятся различные траектории. В Mathcad можно реализовать эту задачу, например, в форме алгоритма, приведенного в листинге 11.10 для решения системы уравнений автокаталитической химической реакции с диффузией. Эта модель, называемая моделью брюсселятора, предложена в 1968 г. Лефевром и Пригожиным. Неизвестные функции отражают динамику концентрации промежуточных продуктов некоторой реальной химической реакции. Параметр модели в равен исходной концентрации катализатора.

Листании 11.10. Построение фазового портрета для модели брюсселятора

Предложенный алгоритм формирует из отдельных матриц решений системы ОДУ с разными начальными условиями объединенную матрицу и. Пары начальных условий задаются в первой строке листинга в виде матрицы v размера 2х10. Последнее означает построение десяти траекторий. Для того чтобы поменять количество траекторий, измените соответствующим образом размер этой матрицы. Затем (рис. 11.14) элементы матрицы и выводятся на график в виде отдельных точек. Отсутствие соединения точек линиями является недостатком алгоритма, но это плата за возможность представить в Mathcad несложным образом сразу большое количество траекторий на фазовой плоскости.

Рис. 11.14. Фазовый портрет брюсселятора при в=0.5 (листинг 11.10)

Как видно из рис. 11.14, все траектории, вышедшие из разных точек, асимптотически стремятся к одному и тому же аттрактору (1,0.5). Из теории динамических систем нам известно, что такой аттрактор называется узлом (с узлом мы уже встречались в примерах разд. 11.1). Конечно, в общем случае при анализе фазового портрета желательно "прощупать" большее число траекторий, задавая более широкий диапазон начальных условий. Не исключено, что в других областях фазовой плоскости траектории будут сходиться к другим аттракторам

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

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

Рис. 11.15. Фазовый портрет брюсселятора при в=2.5

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


8.5. Жесткие системы ОДУ

До сих пор мы имели дело с "хорошими" уравнениями, которые надежно решались численными методами Рунге-Кутты. Однако имеется класс так называемых жестких (stiff) систем ОДУ, для которых стандартные методы практически неприменимы, поскольку их решение требует исключительно малого значения шага численного метода. Некоторые из специальных алгоритмов, разработанных для этих систем, реализованы в Mathcad.

8.5.1. Что такое жесткие ОДУ?

Строгого общепринятого математического определения жестких ОДУ нет. В рамках этой книги будем считать, что жесткие системы — это те уравнения, решение которых получить намного проще с помощью определенных неявных методов, чем с помощью явных методов (типа тех, что были рассмотрены в предыдущих разделах). Примерно такое определение было предложено в 1950-х годах классиками в этой области Кертиссом и Хирш-фельдером. Начнем разговор о жестких ОДУ с примера нежесткого уравнения (листинг 11.11), решение которого показано на рис. 11.16. На том же графике показано решение подобного ОДУ с другим коэффициентом при правой части, равным не -10, а -30. Решение обоих уравнений не составило труда, и численный метод Рунге-Кутты дал правильный результат.

Листании 11.11. Решение нежесткого ОДУ

Рис. 11.16. Решение нежестких ОДУ методом Рунге-Кутты (листинг 11.11)

На рис. 11.17 показано решение того же ОДУ с коэффициентом -50. Вас, несомненно, должен насторожить результат, выданный Mathcad. Характерная "разболтка" решения говорит о неустойчивости алгоритма. Первое, что можно сделать, — увеличить количество шагов в методе Рунге-Кутты. Для этого достаточно добавить третий параметр step в функцию odesoive(t,i,step). После нескольких экспериментов можно подобрать такое значение step, которое будет обеспечивать устойчивость решения. Читатель может самостоятельно убедиться, что при step>20 "разболтка" пропадает, и решение становится похожим на графики, показанные на рис. 11.16.

Рис. 11.17. Неверное решение более жесткого ОДУ методом Рунге-Кутты

Таким образом, во-первых, мы выяснили, что одни и те же уравнения с разными параметрами могут быть как жесткими, так и нежесткими. Во-вторых, чем жестче уравнение, тем больше шагов в обычных численных методах требуется для его устойчивого решения. С классическим примером ОДУ из листинга 11.11 все получилось хорошо, т. к. оно было не очень жестким, и небольшое увеличение числа шагов разрешило все проблемы. Для решения обычными методами более жестких уравнений требуются миллионы, миллиарды и даже большее число шагов.

Некоторые ученые замечают, что в последние годы методы Рунге-Кутты стали уступать свое главенствующее положение среди алгоритмов решения ОДУ методам, способным решать жесткие задачи.

Исторически, интерес к жестким системам возник в середине XX века при изучении уравнений химической кинетики с одновременным присутствием очень медленно и очень быстро протекающих химических реакций. Тогда неожиданно оказалось, что считавшиеся исключительно надежными методы Рунгй-Кутты стали давать сбой при расчете этих задач. Рассмотрим классическую модель взаимодействия трех веществ (Робертсон, 1966), которая как нельзя лучше передает смысл понятия жесткости ОДУ.

Пусть вещество "0"- медленно превращается в "1": "0" -> "1" (со скоростью 0.1), вещество "1" при каталитическом воздействии самого себя превращается очень быстро в вещество "2" : "1" + "1" -> "2" + "1" (10^3). И, наконец, подобным! образом (но со средней скоростью) реагируют вещества "2" и "1" : "1" + "2" -> "0" + "2" (10^2). Система ОДУ, описывающая динамику концентрации реагентов, с попыткой решения методом Рунге-Кутты, приведена в листинге 11.12.

Листинг 11.12. Жесткая сисема ОДУ химической кинетики

Бросается в глаза сильно различающийся порядок коэффициентов при разных слагаемых. Именно степень этого различия чаще всего и определяет жесткость системы ОДУ. В качестве соответствующей характеристики выбирают матрицу Якоби (якобиан) векторной функции Fit,у), т.е. функциональную матрицу, составленную из производных F(t,y) (см. разд. "Частные производные" гл. 7). Чем вырожденнее матрица Якоби, тем жестче система уравнений. В приведенном примере определитель якобиана и вовсе равен нулю при любых значениях у„, yi и у2 (листинг 11.13, вторая строка). В первой строке листинга 11.13 приведено напоминание способа вычисления якобиана средствами Mathcad на примере определения элементов его первой строки.

Листинг 11.13. Якобиан рассматриваемой системы ОДУ химический киметики

Для примера, приведенного в листинге 11.12, стандартным методом Вунге-Кутгы все-таки удается найти решение (оно показано на рис. 11.18). Однако для этого требуется очень большое число шагов, м=2000о, что делает (расчеты очень медленными. При меньшем числе шагов численному алгоритму не удается найти решение. В процессе работы алгоритма оно расходятся, и Mathcad вместо результата выдает ошибку о превышении предельно большого числа.

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

В принципе, можно было бы снизить жесткость системы "вручную", применяя масштабирование. Для этого нужно искусственно уменьшить искомую функцию yl, к примеру, в тысячу раз, разделив все слагаемые в системе ОДУ, содержащие yi, на 1000. После масштабирования для решения полученной системы методом Рунге-Кутты будет достаточно взять всего М=20 шагов.

Рис. 11.18. Решение жесткой системы ОДУ химической кинетики методом Рунге-Кутты (листинг 11.12)


8.5.2. Функции для решения жестких ОДУ

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

  • Radau(y0,t0,t1,M,F) — алгоритм RADAUS для жестких систем ОДУ;
  • stiffb(y0,t0,t1,M,F,j) — алгоритм Булирша-Штера для жестких систем ОДУ;
  • stiffr(y0,t0,t1,M,F,j) — алгоритм Розенброка для жестких систем ОДУ;
    • у0 — вектор начальных значений в точке to;
    • t0.t1 — начальная и конечная точки расчета;
    • M— число шагов численного метода;
    • F — векторная функция F(t,y) размера IXN, задающая систему ОДУ;
    • J- матричная функция J(t,y) размера (N+DXN, составленная из вектора производных функции F(t,y) no t (правый столбец) и ее якобиана (N левых столбцов).

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

Покажем действие этих алгоритмов на том же примере жесткой системы ОДУ химической кинетики (листинг 11.14). Обратите внимание, как следует представлять в данном случае якобиан, сравнив задание матричной функции в предпоследней строке листинга 11.14 с заданием якобиана из листинга 11.13.

Листинг 11.14. Решение жесткой системы ОДУ химической кинетики

Расчеты показывают, что для получения того же результата (см. рис. 11.18) оказалось достаточно в тысячу раз меньшего количества шагов численного алгоритма, чем для стандартного метода Рунге-Кутты! Примерно во столько же раз требуется меньше компьютерного времени на проведение расчетов. Стоит ли говорить, что, если Вы имеете дело с жесткими (в той или иной степени) системами, применение описанных специальных алгоритмов просто необходимо.

Важно заметить, что до сих пор мы имели дело с примером не очень жесткой системы. Попробуйте вместо скоростей упомянутых химических реакций (см. разд. 11.5.1), 0.1, 103 и 102 взять другие числа, например 0.05, 10 и 107, соответственно. Заметим, что такое соотношение скоростей часто встречается в прикладных задачах химической кинетики и определяет куда более жесткую систему ОДУ. Ее уже никак не удается решить стандартными методами, поскольку число шагов численного метода должно быть просто гигантским. А между тем, алгоритмы для жестких ОДУ справляются с этой задачей с легкостью (рис. 11.19), причем практически при тех же значениях шага, что были взяты в листинге 11.14. Обратите внимание, что порядки величины решений для концентраций различных веществ на рис. 11.19 различаются еще сильнее, чем в предыдущем (менее жестком) примере.

Это еще раз доказывает, что одна и та же система ОДУ с различными коэффициентами может быть жесткой в разной степени В частности, приведенный выше пример генератора Ван дер Поля с параметром ц=5000 — это уже пример жесткой задачи.

В заключение приведем соответствующие встроенные функции, которые применяются для решения жестких систем ОДУ не на всем интервале, а только в одной заданной точке t1.

  • radau(y0, t0, t1,acc,F,k,s) —метод RADAUS
  • stiffb(y0,t0,t1,acc,F,j,k,s) — метод Булирша-Штера
  • stiffr(y0,t0,t1,acc,F,j,k,s) —метод Розенброка

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

Рис. 11.19. Решение более жесткой системы ОДУ химической кинетики методом Розенброка

Modifié le: jeudi 11 août 2011, 12:32