Y' = F(x, Y), Y(x0) = Y0

Y(x0) = Y0 - вектор начальных условий;

Y'=(y'1, y'2, …, y'n) - вектор первых производных;

F(x, Y) = (y2, y3, …, yn, f(x,y1, … , yn) - вектор правых частей;

Y = (y2, y3, …, yn) - вектор искомого решения.

Эта система получается в результате следующей замены:

 ,где

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

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

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

Вычислительный блок для решения одного ОДУ состоит из трех частей:

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

- ОДУ и начальные условия, записанные с помощью логического равенства;

- встроенная функция Odesolve(x, b) относительно независимой переменной x на интервале [a, b]; b - верхняя граница отрезка интегрирования. Допустимо и даже предпочтительнее задание функции Odesolve(a, b, step) с тремя параметрами, где step - внутренний параметр численного метода, определяющий количество шагов; чем больше step, тем с лучшей точностью будет получен результат, но тем больше времени будет затрачено на его поиск.

Функция Odesolve возвращает решение задачи в виде функции. Эта функция не имеет символьного представления и может только вернуть численное значение решения уравнения в любой точке интервала интегрирования.

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

Альтернативный метод решения ОДУ заключается в использовании одной из встроенных функций: rkfixed, Rkadapt, или Bulstoer. Все они решают задачу Коши для системы дифференциальных уравнений первого порядка, но каждая из них использует для этого свой метод. Для простых систем не играет большой роли, какой метод использовать - все равно получите решение достаточно быстро и с высокой точностью. Но для сложных или специфических систем бывает, что некоторые методы вообще не могут дать удовлетворительного решения за приемлемое время. Именно для таких сложных, но не редких случаев в MathCAD и введено несколько различных методов решения систем ДУ.

- rkfixed - метод Рунге-Кутта с фиксированным шагом интегрирования. Самый простой и быстрый метод, но далеко не всегда самый точный. Полностью аналогичен использованию функции Odesolve с выбранным в контекстном меню методом Fixed.

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

- Bulstoer - метод Булирша - Штера. Этот метод более эффективен, чем метод Рунге-Кутта, в случае если решение является плавной функцией.

Имена функций Rkadapt и Bulstoer начинаются с прописной буквы. В MathCAD для некоторых имен функций неважно, с какой буквы они записаны, но для перечисленных функций это принципиально, т.к. в MathCAD также существуют функции с такими же именами, только записанные с маленькой буквы - rkadap, bulstoer. Эти функции используются в тех случаях, когда важным является решение задачи в конечной точке интервала интегрирования.

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

Применение встроенных функций в документах MathCAD выглядит сходным образом, т.е. функции Rkadapt и Bulstoer имеют тот же синтаксис, что и выше приведенная функция rkfixed. Назначение аргументов в этих встроенных функциях следующее:

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

- а - начало отрезка, на котором ищется решение системы (отрезка интегрирования). Именно в этой точке значения неизвестных функций принимаются равными элементам вектора y.

- b - конец отрезка интегрирования.

- n - количество частей, на которые разбивается отрезок [a, b] при решении системы. Чем больше это число, тем точнее получается решение, но расчет занимает больше времени.

- F(x,y) - векторная функция, элементы которой содержат правые части уравнений системы в нормальной форме (когда левые части - первые производные от соответствующих функций, а в правых частях производные отсутствуют). Аргументами этой функции являются вектор y, элементы которого соответствуют различным неизвестным функциям системы, и скалярный аргумент x , соответствующий независимой переменной в системе. В случае одного уравнения функция F может быть скалярной функцией, зависящей от двух скалярных переменных x и y.

Возвращаемым значением всех вышеперечисленных встроенных функций является матрица. Первый столбец этой матрицы - это точки, на которые разбивается отрезок [a, b], а остальные столбцы - это значения функций системы в этих точках. Если в аргументе функции rkfixed было указано количество частей n = 100, то матрица будет содержать 101 строку вместе с начальной.

Решение систем обыкновенных дифференциальных уравнений.

Для численного интегрирования систем ОДУ в MathCAD также имеется выбор - либо использовать вычислительный блок Given/Odesolve, либо встроенные функции rkfixed, Rkadapt и Bulstoer.

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

где Y и Y' - соответствующие неизвестные векторные функции переменной t, а F - вектор правых частей системы уравнений первого порядка. Именно векторное представление используется для ввода системы ОДУ в среде MathCAD.

Если в систему ОДУ входят и уравнения высших порядков, то оно тоже сводится к системе уравнений первого порядка, как было показано выше. При этом количество нулевых условий для вычислительного блока Given/Odesolve, а также размер вектора начальных условий y и размер вектора правых частей F(x,y) для встроенных функций rkfixed, Rkadapt и Bulstoer должны быть равны сумме порядков всех уравнений.

Вначале покажем решение систем ОДУ первого порядка с использованием вычислительного блока Given/Odesolve

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

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

Рассмотрим решение этой же системы ОДУ первого порядка с использованием встроенной функции rkfixed.

Полученное решение полностью соответствует вышеприведенному решению с использованием вычислительного блока Given/Odesolve. Следует отметить, что начальные условия здесь задаются в виде вектора y, а функциям x(t) и y(t) соответствуют элементы этого вектора y1 и y2. Вектор начальных условий y и вектор правых частей F имеют размер равный двум, т.к. система состоит из двух уравнений первого порядка. Для системы ОДУ, состоящей из двух уравнений второго порядка, размер этих векторов будет равен четырем

Решение жестких ОДУ и систем ОДУ .

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

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

- Radau - метод Radaus для жестких систем. Полностью аналогичен использованию функции odesolve с выбранным в контекстном меню методом Stiff.

- Stiffb - метод Булирша-Штера, адаптированный для жестких систем.

- Stiffr - метод Розенброка.

Краевые задачи для ОДУ

Постановка краевых задач для ОДУ отличается от задач Коши, рассмотренных выше, тем, что граничные условия для них ставятся не в одной начальной точке, а на обеих границах расчетного интервала. Если имеется система N ОДУ первого порядка, то часть из N условий может быть поставлена на одной границе интервала, а оставшиеся условия - на противоположной границе. В связи с тем, что условия поставлены не на одной, а на обеих границах интервала, краевые задачи нельзя решить изложенными выше методами, предназначенными для задач Коши. Для решения краевой задачи в MathCAD нет отдельной функции. Однако есть функции, позволяющие превратить краевую задачу в задачу Коши. Эти функции «угадывают» недостающие начальные условия, исходя из того, что решение должно удовлетворять заданным условиям в конечной точке интервала интегрирования. Простейшей из функций, предназначенных для приведения краевой задачи к задаче Коши, является функция sbval. Для того, чтобы решить двухточечную краевую задачу с помощью этой функции, следует выполнить следующие действия:

Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13