Функция rkfixed в mathcad как работает
БлогNot. MathCAD: решаем основные типы дифференциальных уравнений встроенными функциями
MathCAD: решаем основные типы дифференциальных уравнений встроенными функциями
Решать дифференциальные уравнения (далее ДУ) в MathCAD, составляя собственные подпрограммы-функции не всегда удобно и экономично по времени, хотя и полезно на этапе обучения. Опишем в этой заметке способы решения основных типов ДУ с помощью стандартных средств пакета, ограничимся простыми примерами.
1. ДУ с разделяющимися переменными. Общая постановка задачи: y’=f(x,y)=g(x)*h(y) , y(x0)=y0 . То есть, f(x,y) допускает представление в виде произведения функций от x и от y .
Для решения уравнения достаточно задать его правую часть как пользовательскую функцию MathCAD, определить интервал поиска решения [x0,x1] , начальное условие y0 и применить стандартную функцию Odesolve . Покажем этот процесс на примере уравнения y’=2x-y+x 2 , x∈[0,2] , y(0)=0 с известным решением y(x)=x 2 :
Решение дифференциальных уравнений с разделяющимися переменными
Знак «равно» в записи уравнений, конечно же, жирный (панель Boolean или сочетание клавиш Ctrl+=).
Функция Odesolve вернула именно функцию y , её нужно смотреть от аргумента, например, y(1)= .
И ещё 2 особенности:
Зная точное решение, графически сравним с ним найденное решение. Как видно на графике, MathCAD справился с задачей отлично.
Графики полученного и точного решения совпадают
2. Неоднородное ДУ первого порядка. В общем виде такое уравнение можно записать как y’=a(x)*y+b(x) . Оно решается аналитически по формуле, которую можно найти в любой книге по решению обыкновенных ДУ:
Формула для решения неоднородного ДУ первого порядка
Здесь С – константа интегрирования. Остаётся применить формулу к конкретному уравнению (возьмём для примера задачу y’+2xy=x*e -x 2 sin(x) , y(0)=1 ) и оценить её символьно:
Аналитическое решение неоднородного ДУ первого порядка в MathCAD
Здесь мы получаем решение в общем виде. Нижний оператор оценён символьно (см. панель «Символика»), а аргумент t используется, так как x в документе «уже занят» (для корректной работы символьной оценки переменные не должны быть определены заранее).
После подстановки начального условия получим частное решение y(t,Y0) , а для проверки решения будет достаточно подставить полученную функцию в исходное уравнение и упростить его символьной функцией simplify . Полученный результат в нашем случае совпал с заданной в условии правой частью. Также для y(x,Y0) , как и для любой функции, можно построить график на нужном интервале изменения x.
Проверка решения неоднородного ДУ первого порядка и построение графика
3. Неоднородное ДУ второго порядка. В общем виде имеем уравнение y» + p(x)*y’ + g(x)*y = f(x) плюс набор краевых условий, количество которых соответствует порядку задачи, например y(0)=. y'(0)=. или y(0)=. y(1)=.
Возьмём уравнение, которое мы мучили вот здесь, решим его стандартными средствами, сравним с известным точным решением и построим график:
Решение неоднородного ДУ второго порядка в MathCAD
Здесь при вызове Odesolve второй параметр, равный единице — это правая граница интервала, третий параметр, равный 10, задаёт количество интервалов. Точное решение u(t) взяли по ссылке. Как видим, даже на 10 интервалах всё очень хорошо совпадает.
4. Система ДУ. Подход к решению системы ДУ покажем на примере. Пусть задана система дифференциальных уравнений
x’ = a*x — y — (x 2 + y 2 )*x,
y’ = a*y + x — (x 2 + y 2 )*y,
x(0)=0, y(0)=1, a=-0.2
Чтобы решить эту систему стандартной функцией rkfixed , нужно задать для неё вектор начальных значений x = (x0, y0) и вектор правых частей D(t,x) .
После этого задача решится вызовом rkfixed , второй и третий параметры ( 0, 20 ) задают интервал по времени t , на котором ищется решение, четвёртый параметр 100 означает количество точек на интервале.
Функция вернёт матрицу решений системы, в которой количество строк соответствует количеству точек на интервале, а количество столбцов — количеству уравнений в системе.
Для построения графика достаточно отобразить зависимость столбцов Zi,1 , Zi,2 от Zi,0 , i=0..99 :
Решение системы ДУ в MathCAD функцией rkfixed
Иллюстрированный самоучитель по MathCAD 11
В Mathcad 11 имеются три встроенные функции, которые позволяют решать поставленную в форме (2-3) задачу Коши различными численными методами.
- rkfixed(y0, t0, t1, M, D) – метод Рунге-Кутты с фиксированным шагом,
- Rkadapt(y0, t0, t1, M, D) – метод Рунге-Кутты с переменным шагом;
- Buistoer(y0, t0, t1, M, D) – метод Булирша-Штера;
- у0 – вектор начальных значений в точке t0 размера 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. В-третьих, точно такой же размер должен быть и у вектора начальных значений у0 (он определен во второй строке листинга).
Не забывайте, что векторную функцию D(t,y) следует определять через компоненты вектора у с помощью кнопки нижнего индекса (Subscript) с наборной панели Calculator (Калькулятор) или нажатием клавиши [. В третьей строке листинга определено число шагов, на которых рассчитывается решение, а его последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (0, 50).
Как выглядит все решение, показано на рис. 11.3. Размер полученной матрицы будет равен (M+DX(N+I), т.е. Ю1хз. Просмотреть все компоненты матрицы и, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки. Как нетрудно сообразить, на этом рисунке отмечено выделением расчетное значение первого искомого вектора у0 на 12-м шаге u12.1=0.01. Это соответствует, с математической точки зрения, найденному значению У0(6.0)=0.07. Для вывода элементов решения в последней точке интервала используйте выражения типа Uм1=7.523×103.
Рис. 11.3. Матрица решений системы уравнений (листинг 11.4)
MathCAD — это просто! Часть 22. Другие методы численного решения дифференциальных уравнений
В прошлый раз мы с вами обсуждали вопросы решения дифференциальных уравнений в мощнейшей из математических сред (то есть в MathCAD) с помощью встроенной функции Odesolve. Функция эта, как мы с вами успели выяснить, в общем-то, всем хороша, однако, как я уже успел в прошлый раз обмолвиться, она предлагает не единственный из всех возможных путей численного решения дифференциальных уравнений в MathCAD. Для того, чтобы чувствовать себя в MathCAD’е как рыба в воде, нужно в совершенстве овладеть всеми способами решения дифференциальных уравнений. Ну, а для начала, наверное, стоит познакомиться с этими самыми способами с помощью вот этой статьи. Кроме того, мы с вами разберемся с такими вопросами, как решение систем «диффуров», поскольку системы дифференциальных уравнений встречаются едва ли не чаще, чем уравнения одиночные. В общем, планы, как видите, прямо-таки наполеоновские, а потому не буду затягивать вступления, а перейду к самому рассказу о работе с дифференциальными уравнениями в MathCAD’е.
Выбор метода интегрирования
Если вы еще помните наш с вами разговор о численном интегрировании в среде MathCAD, то должны помнить и о том, что можно выбирать разные методы интегрирования в зависимости от задачи. В общем-то, поскольку решение дифференциальных уравнений — это, по большому счету, то же самое интегрирование, разумно было бы ожидать, что и при решении ДУ можно будет выбирать метод интегрирования. И действительно, MathCAD (как, впрочем, и всегда) не обманывает наших ожиданий: выбирать способ вычисления решения дифференциального уравнения при использовании функции Odesolve и можно, и нужно.
Хитрости в изменении метода интегрирования для Odesolve никакой нет. Достаточно кликнуть правой кнопкой мыши по названию функции и выбрать в появившемся меню один из трех вариантов: Fixed, Adaptive или Stiff. Думаю, имеет смысл все же сделать некоторые пояснения по поводу того, чем различаются данные способы интегрирования дифференциальных уравнений. Fixed использует фиксированный шаг интегрирования, а Adaptive, как можно понять уже из его названия, — адаптивный (то есть изменяющийся). Если вы не очень хорошо представляете себе, что это все означает, и как они друг от друга отличаются, то имеет смысл вернуться к той статье про численное интегрирование в MathCAD’е, в которой я довольно подробно об этом всем рассказывал. Что касается метода Stiff, то он довольно специфичен и используется для так называемых жестких систем. Если вы работаете с ними, то вы сами знаете, что они жесткие, и вам следует устанавливать флажок именно напротив этого пункта меню. Если же, скажем так, в жесткости системы возникают сомнения, то лучше всего использовать метод Adaptive, потому что он гарантирует более высокую точность решения и избавляет пользователя от хлопот, связанных с подбором правильного шага интегрирования.
Использование функции rkfixed
Как я уже говорил во вступлении к статье, функция Odesolve является не единственным возможным вариантом применяемой для решения дифференциальных уравнений функции. Сейчас мы с вами рассмотрим функцию rkfixed, которая также может быть успешно использована для их решения. Вы можете спросить, зачем нужна еще одна функция для численного решения дифференциальных уравнений, если Odesolve и так самым что ни на есть замечательным образом справляется со своими прямыми обязанностями. Что ж, вопрос резонный, но на самом деле хронологически Odesolve появилась в MathCAD’е позже, чем некоторые другие функции для решения дифференциальных уравнений (в том числе rkfixed). В общем-то, пожалуй, можно было бы и не говорить о более старых функциях, поскольку Odesolve более проста для пользователя в силу своей универсальности, однако бывают случаи, когда приходится сталкиваться и со старыми функциями, которые в новых версиях MathCAD оставлены, в основном, из соображений обратной совместимости со старыми версиями этой математической среды. Поэтому, если вам придется столкнуться с документами, сделанными в старых версиях MathCAD’а, в которых используется для решения дифференциальных уравнений функция rkfixed, то вы не будете пребывать в растерянности, а сможете разобраться, что именно хотел вычислить автор данного документа.
Итак, rkfixed. Пользоваться этой функцией, по сути дела, нужно точно так же, как и Odesolve. То есть сначала должны идти начальные условия для уравнения, а также само уравнение, оформленное в виде функции двух переменных (т.е. dy/dx = f(x, y)). К сожалению, функция rkfixed позволяет решать уравнения исключительно такого вида, а потому, например, работать с уравнением второго порядка вам придется с ее помощью уже как с системой двух дифференциальных уравнений. Помимо этого, нужно указать количество шагов интегрирования, для которых будут рассчитываться значения решения нашего дифференциального уравнения. Использования функции rkfixed вы можете на соответствующей иллюстрации к статье. Из нее, кстати, несложно понять, какие аргументы мы передаем в функцию для решения ДУ. Первый параметр — это начальное значение аргумента, по которому берется производная; второй и третий — это, соответственно, начало и конец того отрезка, на котором будет производиться интегрирование; четвертый параметр — это количество шагов интегрирования, ну, а последний — собственно сама функция f(x, y) для нашего уравнения dy/dx = f(x, y).
Как видите, использовать Odesolve действительно проще и удобнее. Помимо rkfixed, в MathCAD’е есть и другие функции, которые позволяют решать ДУ не с использованием фиксированного шага, а с адаптивным методом интегрирования или со специальными методами решения некоторых видов уравнений. Найти их в справке MathCAD’а не представит никакой сложности, кроме того, мы еще поговорим немного ниже о них. Однако их, как и rkfixed, рекомендуется по возможности не применять, заменяя на Odesolve.
Решение систем дифференциальных уравнений
Что ж, с одиночными дифференциальными уравнениями мы с вами, в общем-то, можно считать, практически разобрались. Конечно, было бы интересно узнать, как работают алгоритмы, лежащие в основе методов численного решения ДУ в MathCAD’е. Возможно, в следующий раз нам с вами и удастся об этом поговорить. Пока же на очереди вопрос более насущный и важный — решение систем дифференциальных уравнений. Поспешу вас обрадовать: мы с вами уже почти умеем решать системы дифференциальных уравнений — недаром же мы с вами останавливались на вопросах использования такой, казалось бы, устаревшей и не сильно для нас с вами полезной функции, как rkfixed. Оказывается, именно она поможет нам решать системы дифференциальных уравнений, с которыми, к сожалению, функция Odesolve никаких дел иметь не хочет. Использование функции rkfixed (и некоторых других) для решения систем ДУ имеет самые что ни на есть минимальные отличия от ее же использования в случае одиночных уравнений. По сути дела, для самого MathCAD’а отличий и вовсе никаких нет — просто скалярные переменные в этой функции заменяются на векторные, с которыми она работает точно так же, как и со скалярными. При этом, правда, действует все то же ограничение на вид уравнений, входящих в систему: они должны быть линейными и не содержать производных выше первой. Впрочем, как я уже говорил, ограничение, связанное с производными высших порядков, довольно просто обходится.
Вы легко обнаружите на иллюстрации к статье, демонстрирующей решение системы уравнений с помощью rkfixed, что совсем не так уж много отличий между решением системы и решением одиночного дифференциального уравнения. Кстати, если внимательно присмотреться к решенной на иллюстрации системе, то можно обнаружить, что она как раз и относится к тому «хитрому» способу решения дифференциальных уравнений второго порядка, к которому приходится прибегать при использовании для решения функции rkfixed.
Теперь, пожалуй, стоит рассказать немного о других функциях, аналогичных по сути rkfixed, но работающих несколько иначе. Функция rkfixed (как видно из ее названия) использует для интегрирования фиксированный шаг, что, с одной стороны, обеспечивает приемлемую точность интегрирования, а с другой — его высокую скорость. Но иногда этой точности бывает недостаточно, и тогда имеет смысл применить функцию rkadapt. В отличие от rkfixed, она имеет изменяющийся соответственно скорости изменения интегрируемой функции шаг интегрирования, а потому в случае быстро изменяющихся функций имеет смысл применять именно ее. Впрочем, ее лучше привыкнуть применять всегда, поскольку адаптивный шаг интегрирования позволяет добиться большей точности результатов вычисления решения дифференциального уравнения. Еще одна функция, которая может применяться в процессе решения — Bulstoer. Она позволяет достигать еще большей точности интегрирования, однако применять ее можно только тогда, когда интегрируемая функция будет достаточно гладкой и достаточно медленно изменяющейся. Насколько именно гладкой и медленно изменяющейся должна быть функция, можно попробовать установить экспериментально, сравнив результаты решения с применением каждой из этих трех функций. Как вы понимаете, все они имеют одинаковые параметры, и, поскольку мы уже разобрали параметры rkfixed, то для двух остальных функций мы их разбирать не будем.
Подведем итоги
Что ж, сказано всего в этот раз о численном решении дифференциальных уравнений и их систем было, как всегда, немало, и настало время по традиции подвести итоги. Как видите, решать системы дифференциальных уравнений в MathCAD уже не так просто, как одиночные уравнения, но, в принципе, это все было вполне ожидаемо, и ничего такого особенно трагичного в этом, как мне кажется, нет. Если не бояться их решать, то решать совсем не так уж и сложно — в любом случае, куда как проще, чем с помощью ручки и бумаги. Впрочем, как видите, решение дифференциальных уравнений и их систем сопряжено с рядом довольно-таки тонких моментов, например, касательно того, какой метод интегрирования выбрать. Для того, чтобы не допустить досадной ошибки, отрицательно влияющей на результат, мало просто быть знакомым с разными методами интегрирования — надо еще хорошо представлять себе особенности решаемой с помощью MathCAD задачи и применяемых в ней функций. То есть, несмотря на всю свою мощь, MathCAD все-таки не может освободить вас от знаний высшей математики, потому что это не компьютерный решатель задач, а всего лишь мощный «навороченный» калькулятор. И успехов при работе с MathCAD можно достичь только если подходить к нему именно как к помощнику, а не как к универсальному математическому решателю.
SF, spaceflyer@tut.by
Компьютерная газета. Статья была опубликована в номере 36 за 2008 год в рубрике soft
Функция rkfixed в mathcad как работает
11. Решение обыкновенных дифференциальных уравнений.
Функция для решения дифференциальных уравнений имеет вид:
Z : = rkfixed ( y , x 1, x 2, npoints , D ), где у – вектор начальных значений решений
X 1 – начало отрезка интегрирования
X 2 – конец отрезка интегрирования
Npoints – число точек интегрирования
D – вектор правых частей
Число строк вектора D равно порядку уравнения. При этом дифференциальное уравнение записывается как бы в виде системы. На рис. 18 y 1 — значение первой производной. Следовательно, первая строка вектора D – это первая производная, вторая строка – вторая производная.
Результаты представляются в виде таблицы, где первый столбец число значений независимой переменной, которое равно npoints , второй – значения искомой функции, третий – значения первой производной.
Для построения графика найденной функции необходимо по оси абсцисс отложить первый столбец таблицы, по оси ординат – второй. Столбец указывается верхним индексом (горячие клавиши Ctrl + 6) (см. рис. 18).