Решение дифференциальных уравнений в среде MathCAD

В самом простом виде задача на собственные значения матрицы формулируется следующим образом: требуется найти такие значения , чтобы матричное уравнение имело решение. В таком случае число называют собственным числом матрицы А, а n- компонентный вектор Х, приводящий уравнение с заданным в тождество - собственным вектором. В вышеприведенном примере собственные вектора матрицы А получены в матрице MS. Проверка проведена для первого столбца матрицы MS и соответствующего ему собственного числа 0=5.439.

Решение систем линейных алгебраических уравнений. Этот вопрос является центральным в вычислительной линейной алгебре.

В математике рассматриваются системы линейных уравнений двух видов - однородные и неоднородные.

Неоднородная система уравнений в матричном виде записывается следующим образом: . Здесь А - матрица коэффициентов системы, В - вектор свободных членов, Х - вектор неизвестных системы.

Неоднородная система имеет одно единственное решение, если определитель матрицы отличен от нуля. Для нахождения точного решения неоднородных систем линейных уравнений в линейной алгебре используются три основных метода:

- метод обратной матрицы, он вам уже известен;

- метод исключений Гаусса;

- метод Крамера.

Неоднородная система линейных уравнений в случае равенства ее определителя нулю имеет множество решений, если ранг матрицы системы равен рангу расширенной матрицы системы, либо не имеет решения, если это условие не выполняется. Решить такие системы в MatCADe можно методом Гаусса.

В выше приведенном примере получили систему из трех уравнений с пятью неизвестными, поэтому решение системы будет иметь два свободных параметра (x4, x5).

Однородная система линейных алгебраических уравнений может быть представлена в виде , т.е. правая часть уравнения представляет вектор из нулевых элементов. Как известно, для того чтобы однородная система линейных уравнений имела решение, определитель соответствующей матрицы должен равняться нулю. Это означает, что количество независимых уравнений в системе (т.е. ранг матрицы) меньше, чем количество неизвестных (т.е. порядок матрицы): rank(A) < n. Но вначале нужно выделить в системе эти самые независимые уравнения. Это делается с помощью функции rref, которая с помощью метода исключений Гаусса приводит матрицу к ступенчатому виду.

Дифференциальные уравнения являются основой огромного количества расчетных задач из самых различных областей науки и техники.

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

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

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

- краевые задачи, для которых заданы определенные соотношения сразу на обеих границах интервала. Из дифференциальных уравнений в частных производных есть возможность решать только уравнения с двумя независимыми переменными: одномерные параболические и гиперболические уравнения, такие как уравнения теплопроводности, диффузии, волновые уравнения, а также двухмерные эллиптические уравнения (уравнения Пуассона и Лапласа).

В MathCAD нет универсальной функции для решения дифференциальных уравнений, а есть около двадцати функций для различных видов уравнений, дополнительных условий и методов решения. Эти функции можно найти в библиотеке Insert/Function, категория “Differential Equation Solving (решение дифференциальных уравнений).

Решение Обыкновенных Дифференциальных Уравнений (ОДУ)

F(x,y,y')=0

F - известная функция трех переменных;

x - независимая переменная на интервале интегрирования[a,b];

y - неизвестная функция;

y' - ее производная.

Функция y(x) является решением дифференциального уравнения, если она при всех x[a,b] удовлетворяет уравнению

F(x,y(x),y'(x))=0

График решения y(x) называется интегральной кривой дифференциального уравнения. Если не заданы начальные условия, таких решений y(x) будет множество. При известных начальных условиях y(x0)= y0 решение y(x) будет единственным. Вычислительный процессор MathCAD может работать только с нормальной формой ОДУ. Нормальная форма ОДУ - это ОДУ, разрешенное относительно производной y'=f(x,y)

Обыкновенным дифференциальным уравнением n-го порядка называется уравнение вида

F(x,y,y',y'', …,y⁽ⁿ⁾)=0

F - известная функция n+2 переменных;

x - независимая переменная на интервале интегрирования[a,b];

y - неизвестная функция;

n - порядок уравнения.

Функция y(x) является решением дифференциального уравнения, если она при всех x[a,b] удовлетворяет уравнению

F(x, y(x), y'(x), y''(x),…, y⁽ⁿ⁾(x))=0

Нормальная форма ОДУ высшего порядка имеет вид

Y⁽ⁿ⁾ =f(x, y, y', …, y^(n-1))

Если не заданы начальные условия, то дифференциальное уравнение n - го порядка имеет бесконечное множество решений, при задании начальных условий y(x0)= y0, y'(x0)= y0,1, y''(x0)= y0,2, …, y^(n-1)(x0)= y0,n-1 решение становится единственным (задача Коши).

Задача Коши для дифференциального уравнения n - го порядка может быть сведена к задаче Коши для нормальной системы n дифференциальных уравнений 1 го порядка, которая в векторной форме имеет вид

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 имеют размер равный двум, т.к. система состоит из двух уравнений первого порядка. Для системы ОДУ, состоящей из двух уравнений второго порядка, размер этих векторов будет равен четырем

1. Поясните работу команд панели Matrix - скалярное и векторное произведение, детерминант матрицы, сумма элементов вектора, операция векторизации.

2. Перечислите три основные группы матричных функций. Расскажите о матричных функциях, возвращающих числовые характеристики. Приведите примеры.

3. Матричные функции, реализующие генерацию матриц и операции работы с блоками матриц.

4. Перечислите матричные функции, реализующие численные алгоритмы решения задач линейной алгебры. Объясните, как работают функции rref и rank.

5. Какие функции вычисляют собственные вектора и собственные числа квадратной матрицы?

6. Решение в системе MathCAD неоднородных систем линейных уравнений, когда определитель матрицы не равен нулю. Три способа.

7. Как осуществляется в системе MathCAD решение неоднородных систем линейных уравнений, когда определитель равен нулю и при условии, что ранг матрицы системы равен рангу расширенной матрицы системы?

8. Как осуществляется в системе MathCAD решение однородных систем линейных уравнений, когда определитель матрицы равен нулю (т.е. ранг матрицы должен быть меньше порядка матрицы)?

9. Какие дифференциальные уравнения называются ОДУ первого порядка? Высшего порядка? Что такое нормальная форма ОДУ первого и высшего порядка? К чему сводятся ОДУ высшего порядка при решении?

10. Можно ли решить дифференциальные уравнения в MathCADе символьно?

11. Как решаются ОДУ с помощью вычислительного блока Given/Odesolve? Какой метод решения реализует функция Odesolve? Как можно изменить метод решения для этой функции?

12. Как решаются ОДУ с помощью встроенной функции rkfixed? Чем функция rkfixed отличается от функции Rkadapt?

13. Как осуществляется решение системы ОДУ с помощью вычислительного блока Given/Odesolve? Приведите примеры.

14. Как осуществляется решение системы ОДУ с помощью функции rkfixed? Приведите примеры.