Решение дифференциальных уравнений в среде MathCAD
p align="left">Перемещать линии ввода внутри формульного блока можно также с помощью клавиш со стрелками или щелкая в нужном месте мышью.Большинство операций редактирования формул реализованы естественным образом, однако некоторые из них несколько отличаются от общепринятых, что связано с особенностью MathCAD как вычислительной системы, например, вставка операторов. Операторы могут быть унарными (действующими на один операнд)и бинарными (действующими на два операнда) . При вставке нового оператора в документ MathCAD определяет, сколько операндов ему требуется. Если в точке вставки оператора один или оба операнда отсутствуют, MathCAD автоматически помещает рядом с оператором один или два местозаполнителя. То выражение в формуле, которое выделено линиями ввода в момент вставки оператора, становится его первым операндом
. Как видно, MathCAD сам расставляет, если необходимо, скобки, чтобы часть формулы, отмеченная линиями ввода, стала первым операндом. Местозаполнители (черный квадратик - для символа и прямоугольная рамка - для оператора) появляются внутри незавершенных формул . Символ в черный квадратик вводится обычным образом, а чтобы в прямоугольную рамку ввести оператор, например +, необходимо курсор расположить перед этой прямоугольной рамкой . Вопросы 1. Назначение системы MathCAD. Какие еще пакеты математических программ вы знаете? 2. Интерфейс пользователя в системе MathCAD. 3. Документ в системе MathCAD(заголовок, расширение при сохранении на диск, типы и расположение блоков, точка привязки блока, размеры блоков, сквозная передача данных в документе). 4. Объясните, что такое входной язык системы MathCAD и язык реализации системы MathCAD? 5. Перечислите основные объекты входного языка системы MathCAD. Расскажите об алфавите языка и о встроенных и пользовательских функциях системы MathCAD. Что такое определение функции и обращение к функции? 6. Константы и переменные в системе MathCAD? Как задаются типы данных в MathCAD? Что такое глобальное и локальное присваивание переменных в документе MathCAD? Как вставляется мнимая единица для комплексных чисел? Что такое ранжированная переменная и как она задается? 7. Как задаются массивы в MathCADе? Как можно добавлять строки и столбцы в готовые матрицы? Как удаляются строки и столбцы из матриц? 8. Перечислите операторы входного языка системы MathCAD? 9. Как осуществляется вывод результатов в системе MathCAD? Как можно настроить формат вывода результатов? Как осуществляется управление процессом вычислений в системе MathCAD? 10. Как работать с единицами измерений физических величин в системе MathCAD? 11. Подробно охарактеризуйте текстовые, графические и математические блоки. Лекция №2. Задачи линейной алгебры и решение дифференциальных уравнений в среде MathCAD В задачах линейной алгебры практически всегда возникает необходимость выполнять различные операции с матрицами. Панель операторов с матрицами находится на панели Math.
Операторы , вам уже знакомы. Напомним только, что оператор вычисляет только детерминант матрицы, а модуль вектора, который равен квадратному корню из суммы квадратов его элементов, вычисляется с помощью оператора , который расположен на панели Calculator. К сожалению, по внешнему виду они не отличаются.
При попытке вычислить модуль вектора с панели Matrix будет ошибочное состояние. Точно также будет ошибочное состояние при попытке вычислить детерминант матрицы с панели Calculator. Рассмотрим неизвестные вам до сих пор операторы панели Matrix. Скалярное произведение векторов определяется как скаляр, равный сумме попарных произведений соответствующих элементов (идентичен обычному оператору умножения). Векторы должны иметь одинаковый размер. Для обозначения скалярного произведения используется символ «точка». Векторное произведение двух векторов u и v с углом между ними равно вектору с модулем , направленным перпендикулярно плоскости векторов u и v. Векторное произведение векторов применимо только для трехкомпонентных векторов. Обозначают векторное произведение символом х, который можно ввести нажатием кнопки на панели Matrix/ - сумма элементов вектора.
- оператор векторизации. Он позволяет провести однотипную операцию над всеми элементами массива (т.е. матрицы или вектора), упрощая тем самым программирование циклов. Например, иногда требуется умножить каждый элемент одного вектора на соответствующий элемент другого вектора. Непосредственно такой операции в MathCAD нет, но ее легко осуществить с помощью векторизации. Оператор векторизации можно использовать только с векторами и матрицами одинакового размера.
Для решения задач линейной алгебры в MathCAD встроены матричные функции. Их можно разделить на три основные группы: - функции определения (генерации) матриц и операции с блоками матриц; - функции вычисления различных числовых характеристик матриц; - функции, реализующие численные алгоритмы решения задач линейной алгебры. Из каждой группы приведем по несколько, наиболее часто используемых функций. Первая группа: 1. matrix(m, n, f) - создает и заполняет матрицу размерности m x n, элемент которой, расположенный в i-ой строке и j-ом столбце равен значению f(i, j) функции f(x, y); 2. diag(v) - создает диагональную матрицу, элементы главной диагонали которой хранятся в векторе v; 3. identity(n) - создает единичную матрицу порядка n; 4. augment(A, B) -объединяет матрицы A и B; матрица B располагается справа от матрицы A, при этом матрицы должны иметь одинаковое число строк; 5. stack(A, B) - объединяет матрицы A и B, матрица В располагается внизу под матрицей А, при этом матрицы должны иметь одинаковое число столбцов; 6. submatrix(A, ir, jr, ic, jc) - формирует матрицу, которая является блоком матрицы А, расположенным в строках с ir по jr и в столбцах с ic по jc, причем ir jr, ic jc.
Вторая группа: 1. last(v) - вычисляет номер последнего элемента вектора V; 2. length(v) - вычисляет количество элементов вектора V; 3. min(v), max(v) - вычисляет минимальное и максимальное значения вектора V; 4. Re(v) - создает вектор из реальных частей комплексных элементов вектора V; 5. Im(v) - создает вектор из мнимых частей комплексных элементов вектора V; 6. sort(V) - сортировка элементов вектора V по возрастанию; 7. reverse (sort(v)) - сортировка элементов вектора V по убыванию; 8. csort (A,n) - сортировка элементов n - го столбца матрицы А по возрастанию (перестановкой строк); 9. rsort (A,n) - сортировка элементов n - ой строки матрица А по возрастанию (перестановкой столбцов); 10. rows(A) - вычисляет число строк в матрице А; 11. cols(A) - вычисляет число столбцов в матрице А; 12. max(A), min(A) - определяет максимальное и минимальное значения матрицы А; 13. tr(A) - вычисляет след квадратной матрицы А (след матрицы равен сумме ее диагональных элементов по главной диагонали); 14. mean(A) - среднее значение элементов матрица А. Действие функций второй группы ясно из их названия, поэтому примеры для них приводить не будем. Третья группа: 1. rref(A) - приведение матрицы к ступенчатому виду с единичным базисным минором (выполняются элементарные операции со строками матрицы: перестановка строк, умножение строки на число, сложение строк); 2. rank(A) - вычисляет ранг матрицы А (количество линейно-независимых строк или это число ненулевых строк ступенчатой матрицы rref(A)); 3. eigenvals(A) - вычисление собственных значений квадратной матрицы А; 4. eigenvecs (A) - вычисление собственных векторов квадратной матрицы А, значением функции является матрица, столбцы которой есть собственные векторы матрицы А, причем порядок следования векторов отвечает порядку следования собственных значений, вычисленных с помощью функции eigenvals(A); 5. eigenvec(A,e) - вычисление собственного вектора матрицы А, отвечающего собственному значению e; 6. normi(A) - max - норма, или - норма (infinity norm). в линейной алгебре используются различные матричные нормы, которые ставят в соответствие матрице некоторую скалярную числовую характеристику; 7. lsolve (A,b) - решение системы линейных алгебраических уравнений вида . Функции третьей группы реализуют, как правило, довольно сложные вычислительные алгоритмы. Приведем примеры на использование функций rref и функций для вычисления собственных значений и собственных векторов матрицы. Задача поиска собственных значений и собственных векторов матрицы очень часто встречается в вычислительной практике.
В самом простом виде задача на собственные значения матрицы формулируется следующим образом: требуется найти такие значения , чтобы матричное уравнение имело решение. В таком случае число называют собственным числом матрицы А, а 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(n))=0
F - известная функция n+2 переменных; x - независимая переменная на интервале интегрирования[a,b]; y - неизвестная функция; n - порядок уравнения. Функция y(x) является решением дифференциального уравнения, если она при всех x[a,b] удовлетворяет уравнению
F(x, y(x), y'(x), y''(x),…, y(n)(x))=0 Нормальная форма ОДУ высшего порядка имеет вид
Y(n) =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? Приведите примеры.
Страницы: 1, 2
|