Как решать дифференциальные уравнения в маткад

Как решать дифференциальные уравнения в маткад


/ КР_студентам_cбросить / Пособие_MathCAD / Лаб_5_Решение диф_уравнений

5. Решение дифференциальных уравнений с помощью

встроенных функций MathCAD

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

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

Для численного интегрирования одного ОДУ у пользователя Mathcad имеется выбор - либо использовать вычислительный блок Given – Odesolve( ), либо встроенные функции. Первый путь предпочтительнее из соображений наглядности представления задачи и результатов, а второй дает пользователю больше рычагов воздействия на параметры численного метода. Рассмотрим последовательно оба варианта решения.

Вычислительный блокGiven – Odesolve( )

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

- Given- ключевое слово;

- ОДУ и начальное условие, записанное с помощью Булевых операторов, причем начальное условие должно быть в форме у(t0) = b;

- Odesolve(t, t1)- встроенная функция для решения ОДУ относительно переменнойtна интервале(t0,t1).

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

Пример решения задачи Коши для ОДУ первого порядка у'=у-у 2 посредством вычислительного блока Given – Odesolve( ) приведен на рис. 5.1. Вставлять логические операторы нужно при помощи панели инструментов Boolean (Булевы операторы). При вводе с клавиатуры логического знака равенства нужно использовать сочетание клавиш Ctrl =. Символ производной можно ввести как средствами панели Calculus (Вычисления), как это сделано на рис. 5.1, так и в виде штриха (‘), набрав его с помощью сочетания клавиш Ctrl +F7.

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

Mathcad требует, чтобы конечная точка интегрирования ОДУ лежала правее начальной (t0<t1),иначе будет выдано сообщение об ошибке. ФункцияOdesolve( )возвращает функциюy(t),определенную на интервале(t0,t1).

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

Использование встроенных функцийrkfixed( ), Rkadapt( ), Bulstoer( )

Альтернативным методом решения ОДУ является использование встроенных функций rkfixed( ), Rkadapt( )илиBulstoer( ).Этот метод несколько проигрывает первому и в простоте, и в наглядности, однако иногда приходится решать ОДУ с помощью альтернативного метода, например, при следующих обстоятельствах:

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

- ответ предпочтительнее получить в виде матрицы, а не функции;

- в старых версиях Mathcad’aне было вычислительного блока, а у пользователя имеется много документов, созданных с помощью альтернативного метода.

Приведем пример использования функции rkfixed( ) для решения той же самой задачи Коши для ОДУ первого порядка у' =у-у 2. Здесь необходимо явно задать начальное значение y:=0.01, описать первую производную D(t,y):=y-y 2 и указать количество точек интегрирования m=100. Результат получаем в виде матрицы y размерностью (m+1) x 2. В первом столбце содержатся значения аргумента t в интервале от 0 до 10, во втором столбце - значения искомой функции y(t), рис. 5.2.

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

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

Рассмотренный пример взят из области математической экологии и описывает динамику популяций с внутривидовой конкуренцией. Сначала происходит рост численности популяции, близкий к экспоненциальному, а затем осуществляется выход на стационарное состояние.

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

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

Внутри вычислительного блока Given-Odesolve( ):

- ОДУ должно быть линейно относительно старшей производной;

- начальные условия должны иметь форму y(t0)=bилиy ( n ) (t0)=b;

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

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

Рис. 5.3. Решение дифференциального уравнения второго порядка

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

- rkfixed(y, t0, t1, m, D)- метод Рунге-Кутта с фиксированным шагом;

- Rkadapt(y, t0, t1, m, D)- метод Рунге-Кутта с переменным шагом;

- Bulstoer(y, t0, t1, m, D) - метод Bulirsch-Stoer;

- у - вектор начальных значений в точке t0 размерностью n x1;

- t0- начальная точка интервала;

- t1- конечная точка интервала;

- m- число шагов, на которых численный метод находит решение;

- D - векторная функция дифференциальных уравнений размером n x1 от двух аргументов - скалярного t, по которому взяты производные и векторного у, элементы которого являются младшими производными искомой функции.

Покажем на примере дифференциального уравнения третьего порядка

методику формирования векторной функции D(t,y).Для этой цели нужно разрешить уравнение относительно старшей производной

В результате замены получим:

После этого можно записать векторную функцию D(t,y):

Векторная функцияD(t,y) имеет количество компонентов, равное порядку дифференциального уравнения.

Поскольку младших производных нет, то вместо уравнений записываются их обозначения (y1 иy2 ). Вместо старшей производной записывается правая частьy3. Начальные условия должны задаваться для искомой функции и каждой младшей производной, например:

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

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

Рис. 4.4. Решение системы ОДУ второго порядка

Система ОДУ на рис. 5.4 задается с помощью векторной функцииD(t,y),составленной по описанной выше методике. Такой же размер имеет вектор начальных значенийу. Число шагов, на которых определяется решение, задано параметромm=100. Решение системы ОДУ осуществляется на интервале (0, 50). Просмотреть все компоненты матрицы можно с помощью вертикальной полосы прокрутки.

Решения обыкновенных дифференциальных уравнений часто удобнее изображать не в функции от аргумента t, а в фазовом пространстве. Фазовый портрет системы - это кривая на фазовой плоскости, построенная в координатахy’(t)иy(t),рис. 5.5.

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

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

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

Рис. 5.5. Фазовый портрет решения системы ОДУ

Функция Rkadapt( ) может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо, либо существуют участки медленных и быстрых его изменений. Метод Рунге-Кутта с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений - частыми. В результате, для достижения одинаковой точности требуется меньшее число шагов, чем для rkfixed( ). Метод Bulirsch-Stoer часто оказывается более эффективным для поиска гладких решений.

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

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

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

- rkadapt(y,t0,t1,асе, D, k, s)- метод Рунге-Кутта с переменным шагом;

- bulstoer(y, t0, t1, acc, D, k, s) - метод Bulirsch-Stoer;

- y- вектор начальных значений в точке t0;

- t0,t1- начальная и конечная точки интервала;

- асе- погрешность вычисления (чем она меньше, тем с лучшей точностью будет найдено решение; рекомендуется выбирать значения погрешности в районе 0.001);

- D- векторная функция, задающая систему ОДУ;

- k- максимальное число шагов, на которых численный метод будет находить решение;

- s- минимально допустимая величина шага.

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

Пример использования функции bulstoer( )приведен на рис. 5.6. Структура матрицы результатовuтакая же, как и в случае решения системы ОДУ на рис. 5.4. Однако в данном случае нужно значение последней точки интервала. Поскольку выполненное количество шагов неизвестно, т. е. размер матрицы заранее неизвестен, то его необходимо предварительно определить с помощью функцииlength( ) и выбрать последнюю точку. Для этой цели необходимо транспонировать матрицуu.

Чтобы использовать альтернативный численный метод, достаточно на рис. 5.6. заменить функцию buistoer( )наrkadapt( ).Эти функции не предназначены для нахождения решения в промежуточных точках интервала, хотя они и выдают их в матрице-результате.

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

Некоторые примеры решения ОДУ

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

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

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

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

Динамические системы с таким поведением называют негрубыми.

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

Рис. 5.8. График решения (слева) и фазовый портрет (справа)

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

Рис. 5.9. Модель Ван дер Поля (μ=1)

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

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

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

Рис. 5.11. Модель Лоренца

Рис. 5.12. Аттрактор Лоренца

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

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

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

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

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

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

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

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

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

Жесткие системы ДУ

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

Строгого общепринятого математического определения жестких ДУ нет. Поэтому обычно считают, что жесткие системы - это те уравнения, решение которых получить намного проще с помощью определенных неявных методов, чем с помощью явных методов.

Рис. 5.14. Программа модели брюсселятора

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

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

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

- Radau(y0,t0,t1,M,F)- алгоритм RADAUS для жестких систем ДУ;

- stiffb(y0,t0,t1,M,F,J)- алгоритмBulirsch-Stoerдля жестких систем ДУ;

- stiffr(y0,t0,t1,M,F,J)- алгоритм Розенброка для жестких систем ДУ;

- у- вектор начальных значений в точкеt0;

- t0,t1- начальная и конечная точки расчета;

- m - число шагов численного метода;

- F— векторная функцияF(t,y)размера1хn, задающая систему ДУ;

- J- матричная функцияJ(t,y)размера (n + 1)xn, составленная из вектора производных функцииF(t,y) пot(правый столбец) и ее якобиана (nлевых столбцов).

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

Приведем пример решения жесткой системы ДУ (рис. 5.16).

Рис. 5.16. Решение жесткой системы ДУ

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

Рис. 5. 17. График решения жесткой системы

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

- radau(y, t0, t1,acc,F,k,s) - метод RADAUS;

- stiffb(y,t0,t1,acc,F,j,k,s) - метод Bulirsch-Stoer;

- stiffr(y0,t0,t1,acc,F,j,k,s) - метод Розенброка.

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

Лабораторная работа 5

Цель работы. -изучить функции для решения дифференциальных уравнений;

-решить обыкновенные дифференциальные уравнения первого и второго

- построить таблицы значений и двумерные графики искомых функций;



как решать дифференциальные уравнения в маткад:/ КР_студентам_cбросить / Пособие_MathCAD / Лаб_5_Решение диф_уравнений 5. Решение дифференциальных уравнений с помощью встроенных функций MathCAD Обыкновенное дифференциальное уравнение

как решать дифференциальные уравнения в маткад

Как решать дифференциальные уравнения в маткад 12 8 10