Как решить дифференциальное уравнение в mathcad

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

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

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

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

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

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

— краевые задачи, для которых заданы определенные соотношения сразу на обеих границах интервала.

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

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

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

ОДУ первого порядка .

ОДУ первого порядка называется уравнение

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

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

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

y’ – ее производная.

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

График решения y(x) называется интегральной кривой дифференциального уравнения. Если не заданы начальные условия, таких решений y(x) будет множество. При известных начальных условиях y(x₀)= y₀ решение y(x) будет единственным.

Вычислительный процессор MathCAD может работать только с нормальной формой ОДУ. Нормальная форма ОДУ – это ОДУ, разрешенное относительно производной

ОДУ высших порядков .

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

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

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

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

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

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

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

Y ( n ) =f(x, y, y’, …, y ( n -1) )

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

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

Y(x₀) = Y₀ – вектор начальных условий;

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

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

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

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

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

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

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

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

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

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

1. Задайте вектор v с количеством элементов равным количеству недостающих начальных условий. Значения элементов этого вектора – это начальные приближения, исходя из которых будет происходить поиск недостающих начальных условий. На данном этапе не конкретизируется, какой из элементов вектора будет соответствовать начальному значению той или иной неизвестной функции в задаче.
2. Задайте функцию F (x,y). Эта функция уже описывалась выше. Она представляет собой вектор, каждый элемент которого – это правая часть одного из уравнений системы.
3. Задайте еще одну векторную функцию load (x,v). Это функция от скалярного аргумента x и вектора v, который имеет столько же компонент, сколько недостающих начальных условий в системе. Сам вектор load должен содержать такое же количество элементов, как и вектор F, т.е. столько, сколько должно быть начальных условий в задаче. Если начальное значение какой-либо из функций известно, то соответствующий элемент вектора load должен содержать это значение. Для функций, начальное значение которых неизвестно, соответствующий элемент вектора load должен содержать один из элементов вектора v.
4. Следует задать еще одну некоторую функцию score (x,y). Аргументы этой функции – скаляр x и вектор y, который имеет столько элементов, сколько уравнений в системе. Количество компонент вектора score должно равняться количеству граничных условий, заданных в конечной точке отрезка интегрирования. На самом деле каждая компонента этого вектора задает одно из граничных условий в конечной точке. Например, если в задаче есть граничной условие y_i(b)=c, то один из элементов вектора score должен быть функцией, которая обращается в нуль при значениях x=b и y(b)=c. Конкретный вид этой функции не играет особой роли, поэтому проще всего задавать ее в таком виде: score_k(x,y):=y_i-c. Таким же образом должны быть заданы все элементы вектора score для всех конечных условий задачи.
5. Теперь все введенные величины нужно использовать как аргументы в функции sbval. Использование этой функции выглядит следующим образом: Y:=sbval(v,a,b,F,load,score). Аргументы a и b – это начало и конец отрезка интегрирования.
6. Результатом функции sbval будет вектор, содержащий недостающие начальные значения. Их последовательность задается той последовательностью, в которой были использованы компоненты вектора v в функции load.
7. Постройте вектор начальных значений, используя известные начальные значения, а также элементы вектора Y, там где значения были неизвестны.
8. Теперь можно решить полученную задачу как задачу Коши, с помощью, например, функции Rkadapt.

Ниже приведено решение дифференциального уравнения на отрезке [0,3] с начальным условием y(0)=1 и граничным условием y(3)=4. Недостающее начальное условие – y’(0).

Решение уравнений в частных производных .

Постановка задач для уравнений в частных производных включает определение самого уравнения (или системы нескольких уравнений), а также необходимого количества краевых условий (число и характер которых зависит от специфики уравнения). Согласно своему названию, уравнения должны содержать частные производные неизвестной функции u (или нескольких функций, если уравнений несколько) по различным аргументам, например, пространственной переменной x и времени t. Соответственно, для решения задачи требуется вычислить функцию нескольких переменных, например, u(x,t) в некоторой области определения аргументов 0 £ x £ L и 0£ t £ T. Граничные условия определяются как заданные временные зависимости функции u, или производных этой функции, на границах расчетной области 0 и L, а начальные – как заданная функция u(x,0).

Сами уравнения в частных производных (несколько условно) можно разделить на три основных типа:

— параболические – содержащие первую производную по одной переменной и вторую – по другой, причем все эти производные входят в уравнение с одинаковым знаком;

— гиперболические – содержащие первую производную по одной переменной и вторую – по другой, входящие в уравнения с разными знаками;

— эллиптические – содержащие только вторые производные, причем одного знака.

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

Дифференциальные уравнения параболического типа.

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

1. Ключевое слово given
2. Уравнение, которое нужно решить. Уравнение должно иметь такой вид: . Для ввода производных в данном случае нельзя пользоваться обычным оператором производной, а нужно пользоваться нижним индексом, как это обычно делается в литературе для записи уравнений в частных производных. При этом нижний индекс набирается не как числовой, а как буквенный (точка).
3. Граничные условия для функции u(x,t). Если уравнение второго порядка по x, то граничных условий должно быть два. Можно использовать как граничные условия Дирихле (u(x₀,t) = w(t)), так и граничные условия Неймана (u_x(x₀,t) = w(t)) или их комбинацию, как в вышеприведенном примере.
4. Начальное значение для неизвестной функции – u(x,0).
5. функция pdsolve(u,x,xrange,t,trange,xpts,tpts). Ее аргументы имеют следующее назначение

— u – имя функции, относительно которой решается уравнение. Для системы уравнений здесь должен быть вектор имен функций (как в odesolve).

— X- имя пространственной переменной.

— Xrange – двухкомпонентный вектор, задающий начало и конец интервала изменения пространственной переменной.

— T – имя временной переменно й. Основная разница между пространственной и временной переменными в данном случае – это то, что все уравнения могут содержать только первые производные по временной переменной.

— Trange – еще один двухкомпонентный вектор. Этот вектор задает начало и конец временного интервала, на котором решается задача.

— Xpts, tpts – количество точек, разбивающих для интегрирования пространственный и временной интервалы соответственно. Эти два параметра можно не указывать, тогда количество точек будет выбрано автоматически из соображений достаточной точности. Рекомендуется задавать эти параметры во всех задачах, кроме простых, поскольку во многих случаях высокая точность вычислений теряет смысл из-за погрешности, вносимой самим методом.

Выше приведен пример решения одномерной задачи теплопроводности для одного бруска, один конец которого теплоизолирован, а ругой поддерживается при определенной температуре

Дифференциальные уравнения гиперболического типа.

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

Помимо вычислительного блока given/pdsolve, для решения параболических и гиперболических уравнений можно использовать встроенную функцию numol.

Эллиптические уравнения .

Решение эллиптических уравнений в частных производных реализовано только для единственного типа задач – двумерного уравнения Пуассона. Это уравнение содержит вторые производные функции u(x,y) по двум пространственным переменным. Для его решения в MathCAD есть отдельные встроенные функции multigrid и relax.

Вопросы

1. Какие дифференциальные уравнения называются ОДУ первого порядка? Высшего порядка? Что такое нормальная форма ОДУ первого и высшего порядка? К чему сводятся ОДУ высшего порядка при решении?
2. Можно ли решить дифференциальные уравнения в MathCADе символьно?
3. Как решаются ОДУ с помощью вычислительного блока Given/Odesolve? Какой метод решения реализует функция Odesolve? Как можно изменить метод решения для этой функции?
4. Как решаются ОДУ с помощью встроенной функции rkfixed? Чем функция rkfixed отличается от функции Rkadapt?
5. Как осуществляется решение системы ОДУ с помощью вычислительного блока Given/Odesolve? Приведите примеры.
6. Как осуществляется решение системы ОДУ с помощью функции rkfixed? Приведите примеры.

Понравилась статья? Добавь ее в закладку (CTRL+D) и не забудь поделиться с друзьями:

6.02. Задача Коши для дифференциального уравнения первого порядка

Общее решение любого дифференциального уравнения первого порядка F(х; у; у’) = 0, как это следует из схемы его получения (1.3), содержит бесчисленное множество частных решений. Возникает естественный вопрос: как из этого множество частных решений выделить интересующее нас конкретное частное решение? Иначе говоря, как из множества интегральных кривых данного дифференциального уравнения выделить нужную интегральную кривую?

Ответ почти очевиден: для этого на плоскости Хоу нужно задать некоторую точку , через которую должна пройти искомая интегральная кривая. Тогда её уравнение и будет тем частным решением, которое выделяется из прочих (рис. 6.2).

Задание точки равносильно заданию условия для искомого, выделяемого из прочих, частного решения данного дифференциального уравнения. Это условие называется Начальным условием для дифференциального уравнения первого порядка F(х; у; у’) = 0. Начальным оно называется потому, что очень часто в реальных задачах по исследованию различного рода процессов роль независимой переменной Х играет время T, а начальным значением Х0 является начальный момент времени T0 (обычно T0 = 0). Тогда начальное условие показывает, какое значение У0 имела искомая функция , описывающая исследуемый процесс, в начальный момент времени Х0. Ну, а сама функция , если нас не интересует предистория процесса, то есть времена Х < х0, ищется для Х > х0.

Если дифференциальное уравнение первого порядка F(х; у; у’) = 0 задано вместе с начальным для него условием , То говорят, что для этого уравнения задана Задача Коши:

Решить её — это значит найти те частные решения дифференциального уравнения F(х; у; у’) = 0 , которые еще удовлетворяют и заданному начальному условию . С точки зрения рисунка 6.2 решить задачу Коши (2.1) – это значит найти уравнения всех интегральных кривых дифференциального уравнения F(х; у; у’) = 0, проходящих через начальную точку .

Как правило, задача Коши (2.1) имеет единственное решение . То есть через заданную начальную точку проходит единственная интегральная кривая дифференциального уравнения F(х; у; у’) = 0 (как на рис. 6.2). Но бывает, что задача Коши не имеет решений. То есть бывает, что ни одна из интегральных кривых не проходит через заданную начальную точку . Тогда такая точка называется Особой точкой дифференциального уравнения. А бывает, что задача Коши имеет несколько решений. То есть бывает, что через начальную точку проходит несколько интегральных кривых. Сколько решений будет у задачи Коши (2.1) и каковы они, выясняется в процессе её решения. А Схема решения задачи Коши (2.1) такова:

1. Решаем дифференциальное уравнение F(х; у; у’) = 0 и находим все его решения. То есть находим общее решение (общий интеграл) и возможные особые решения .

2. Подставляем начальные значения Х = х0 И У = у0 в общее решение и находим соответствующее значение (значения) константы С:

3. Подставляем каждое из найденных значений С В общее решение и получаем частные решения

Являющиеся решением задачи Коши. Это те решения этой задачи, которые выделяются из общего решения дифференцированного уравнения F(х; у; у’) = 0.

4. Проверяем, нет ли среди особых решений Дифференциального уравнения F(х; у; у’) = 0 таких, которые удовлетворяют начальному условию У(х0) = у0. Если такие найдутся, они тоже будут решениями задачи Коши (2.1).

Пример1. Решить задачу Коши:

1. Сначала решим дифференциальное уравнение . Оно уже решено ранее – его решение найдено в примере 3, §1:

— общее решение; – особое решение.

2. Подставим начальные значения В общее решение и найдем С:

3. Подставим в общее решение и получим частное решение

Эта функция является решением данной задачи Коши.

4. Обратим внимание на особое решение У=0. Начальному условию У(0)=1 оно не удовлетворяет, поэтому решением данной задачи Коши не является.

Ответ: — единственное решение поставленной задачи Коши.

Пример 2. Материальное тело поднято на высоту H и в начальный момент времени T=0 отпущено в свободное падение. Описать математически процесс падения тела. А именно, найти зависимость ν = ν(T) скорости ν падающего тела от времени T, и найти зависимость S = S(T) пути S, пройденного падающим телом, от времени T. Сопротивлением воздуха пренебречь.

Решение. Как известно, все свободно падающие тела падают с постоянным ускорением G ≈ 9,8 м/сек2 — с ускорением свободного падения. А так как ускорение – это производная от скорости, то получаем: . Это — дифференциальное уравнение первого порядка для искомой функции . Учтём еше, что в начальный момент времени T = 0 тело покоилось, а значит, выполняется начальное условие: . В итоге для определения функции Получаем задачу Коши:

Решим эту задачу.

1. Сначала решим дифференциальное уравнения:

Это – общее решение уравнения , содержащее все его решения. Особых решений у него нет.

2. Используем начальное условие и найдем С:

0 = G0+С => С = 0.

3. Подставим С=0 в общее решение V=Gt+C и получим окончательно: V=Gt. Это и есть решение поставленной задачи Коши (единственное). И заодно V=Gt — это искомая зависимость скорости V падающего тела от времени T.

А теперь займёмся поиском зависимости S=S(t) Пути S От времени T. Учтём, что и что . Тогда для определения этой зависимости получим следующую задачу Коши:

Решим эту задачу.

1. Сначала решим дифференциальное уравнение :

Это — общее решение уравнения , содержащее все его решения. Особых решений у него нет.

2. Используем начальное условие и найдём С:

3. Подставим С=0 в общее решение и получим окончательно: . Это и есть решение рассматриваемой задачи Коши. И заодно — это искомая зависимость пути S, проходимого свободно падающим телом, от времени T.

Ответ: — известные школьные формулы.

Пример3. Дать математическое описание демографического процесса (процесса изменение численности населения со временем) для достаточно крупного населённого региона, если в начальный момент времени численность населения региона составляла человек.

Решение. Пусть – искомая зависимость численности населения региона от времени . И пусть за время , прошедшее с некоторого момента до момента , родилось человек и умерло человек. Эти количества, очевидно, пропорциональны как исходной (в момент ) численности населения , так и величине временного промежутка . То есть

Здесь и – некоторые числовые коэффициенты, связанные соответственно с уровнем рождаемости и уровнем смертности в данном регионе. Тогда общее изменение численности населения за время найдется по формуле:

Здесь . Из получённого равенства следует: . Устремляя здесь (при этом, очевидно, и ), то есть переходя к бесконечно малым и , получим:

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

Решим эту задачу.

1. Сначала решим дифференциальное уравнение . Функция является его очевидным частным решением. Но это, очевидно, не та функция, которую мы ищем – она не удовлетворяет начальному условию, да и вообще она означает, что население в регионе отсутствует.

Будем искать те решения уравнения для которых :

Итак, — общее решение дифференциального уравнения . В него, кстати, при С = 0 входит и отмеченное ранее нулевое решение . То есть в найденном общем решении содержатся все решения дифференциального уравнения.

2. Используем начальное условие и найдём С:

3. Подставим в общее решение и получим искомое решение задачи Коши:

Это и есть искомая зависимость Численности населения региона от времени .

Проанализируем эту зависимость.

а) Если , то численность населения экспоненциально растёт со временем (рис. 6.3(а)).

б) Если , то численность населения Экспоненциально убывает со временем (рис. 6.3(б)).

в) Если , то , то есть численность населения региона не меняется (рис. 6.3(в).

Какой именно будет величина для данного региона, можно выяснить опытным путём. Пусть, например, перепись населения показала, что в некоторый момент времени в регионе проживало человек. Подставляя эти данные в формулу , можем найти :

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

Пример 4. Рассмотрим задачу о математической модели естественного роста выпуска продукции.

Пусть — объем продукции некоторого предприятия, реализованной моменту времени . Будем считать, что вся продукция реализуется по некоторой фиксированной цене за единицу продукции независимо от объема продаж . Это значит, что рынок данной продукции длительное время является насыщенным – удается продавать по фиксированной цене практически любые объемы этой продукции.

Доход от продаж составит: . Будем считать, что некоторая часть этого дохода используется в качестве инвестиций в производство выпускаемой продукции. То есть объем инвестиций составит:

Здесь – так называется Норма инвестиций. Она показывает, какая часть дохода возвращается в производство.

Чем больше объем инвестиций , тем быстрее растёт объем производства . В модели естественного роста это значит, что скорость роста объема производства (так называемая Акселерация производства) пропорциональна объему инвестиций :

— так называемая Норма акселерации, которая показывает, каким должен быть объём инвестиций , чтобы обеспечить единичную скорость роста объема производства (обеспечить рост на единицу продукции за единицу времени). Подставляя (2.4) в (2.5), получим

Где – числовой коэффициент. Равенство (2.7) представляет собой дифференциальное уравнение первого порядка для функции . Дополняя его некоторым начальным условием , получим задачу Коши:

Эта задача полностью совпадает с задачей Коши для демографического процесса (см. пример 3). Значит, у них полностью совпадают и решения:

Заметим, что условие постоянства цены единицы продаваемой продукции, то есть условие насыщенности рынка, не может выполнятся всегда, при любых . С увеличением объема продаж на некотором этапе рынок насыщается, спрос на товар падает, и дальнейшее увеличении объема продаж возможно лишь при снижении цены на него – в соответствии с классической убывающей кривой спроса . Если учесть эту зависимость от , то выражение (2.4) для примет вид:

А вместо (2.7) из (2.5) получим:

Где . Это дифференциальное уравнение вместе с начальным условием составит задачу Коши:

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

1. Сформулировать и решить задачу по определению скорости V=V(T) свободно падающего тела массой M при условии, что учитывается сопротивление воздуха, пропорциональное скорости падения тела.

2. Сформулировать и решить задачу по определению объема У=y(t) реализованной продукции, если известно, что кривая спроса Р= р(у) задаётся уравнением Р=2-у; норма инвестиций M=0,5; норма акселерации ; У(0)=0,5 – начальное условие.

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

• Главная
• Заказать работу
• Стоимость решения
• Варианты оплаты
• Ответы на вопросы (FAQ)
• Отзывы о нас
• Примеры решения задач
• Методички по математике
• Помощь по всем предметам
• Заработок для студентов

Mathcad решаем дифференциальное уравнение 2 порядка

Mathcad решаем дифференциальное уравнение 2 порядка

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

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

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

Дифференциальные уравнения первого порядка

Дифференциальное уравнение первого порядка — это уравнение, которое не содержит производных выше первого порядка от неизвестной функции. На Рисунке 1 показан пример того, как решить относительно простое дифференциальное уравнение:

с начальными условиями: y(0) = 4

Функция rkfixed на Рисунке 1 использует для поиска решения метод Рунге-Кутты четвертого порядка. В результате решения получается матрица, имеющая два следующих столбца:

• Первый столбец содержит точки, в которых ищется решение дифференциального уравнения.
• Второй столбец содержит значения найденного решения в соответствующих точках.

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

Функция rkfixed имеет следующие аргументы:

y =	Вектор начальных условий размерности n, где n — порядок дифференциального уравнения или число уравнений в системе (если решается система уравнений). Для дифференциального уравнения первого порядка, как, например, для уравнения, приведенного на Рисунке 1, вектор начальных значений вырождается в одну точку y0 = y(x1).
x1, x2 =	Граничные точки интервала, на котором ищется решение дифференциальных уравнений. Начальные условия, заданные в векторе y, — это значение решения в точке x1.
npoints =	Число точек (не считая начальной точки), в которых ищется приближенное решение. При помощи этого аргумента определяется число строк (1 + npoints) в матрице, возвращаемой функцией rkfixed.
D (x, y) =	Функция, возвращающая значение в виде вектора из n элементов, содержащих первые производные неизвестных функций.

Наиболее трудная часть решения дифференциального уравнения состоит в определении функции D(x, y), которая содержит вектор первых производных от неизвестных функций. В примере, приведенном на Рисунке 1, было достаточно просто разрешить уравнение относительно первой производной , и определить функцию D(x, y). Иногда, особенно в случае нелинейных дифференциальных уравнений, это может быть трудно. В таких случаях иногда удаётся разрешить уравнение относительно в символьном виде и подставить это решение в определение для функции D(x, y). Используйте для этого команду Решить относительно переменной из меню Символика.

Рисунок 2: Более сложный пример, содержащий нелинейное дифференциальное уравнение.

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

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

• Вектор начальных условий y теперь состоит из двух элементов: значений функции и её первой производной в начальной точке интервала x1.
• Функция D(t, y) является теперь вектором с двумя элементами:

Пример, приведенный на Рисунке 3, показывает, как решить следующее дифференциальное уравнение второго порядка:

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

Уравнения более высокого порядка

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

• Вектор начальных значений y теперь состоит из n элементов, определяющих начальные условия для искомой функции и ее производных y, y’ , y». y (n-1)
• Функция D является теперь вектором, содержащим n элементов:

Пример, приведенный на Рисунке 4, показывает, как решить следующее дифференциальное уравнение четвертого порядка:

с начальными условиями:

Рисунок 4: Решение дифференциального уравнения более высокого порядка.

Исправляем ошибки: Нашли опечатку? Выделите ее мышкой и нажмите Ctrl+Enter

Mathcad решаем дифференциальное уравнение 2 порядка

Электронный курс по MathCAD

5.2 Решение дифференциальных уравнений и систем.(Задача Коши и граничные задачи).

Решение одиночного дифференциального уравнения.

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

Odesolve(x,b,[step]) — Возвращает функцию, которая является решением дифференциального уравнения. Используется в блоке с оператором Given.
x — переменная интегрирования, действительное число
b — конечная точка отрезка интегрирования
step — величина шага по переменной интегрирования (необязательный аргумент)

1. Уравнение должно быть линейным относительно старшей производной.
2. Число заданных начальных или граничных условий внутри блока должно быть равно порядку уравнения.
3. При записи уравнения для обозначения производных функции используйте специальные кнопки с панели Math или ‘ (штрих) — [Ctrl+F7], для знака равенства = [Ctrl+=] (в том числе и для дополнительных условий).
4. Конечная точка должна быть больше начальной.
5. Не допускаются начальные и граничные условия смешанного типа (f ‘(a)+f(a)=5).
6. Искомая функция в блоке дложна быть обязательно с аргументом ( f(x))

Численное решение задачи Коши для дифференциальных уравнений и систем.

Для численного решения задачи Коши для дифференциальных уравнений и систем могут быть использованы функции: rkfixed(y,x1,x2,n,F) — возвращает матрицу решений системы уравнений методом Рунге-Кутта 4-го порядка при фиксированном шаге по x rkadapt(y,x1,x2,n,F) — ищет решение с переменным шагом ( там, где решение меняется медленнее, шаг увеличивается, а в области быстрого изменения решения шаг функции уменьшается). Возвращается решение с равным шагом. Функция работает быстрее, чем rkfixed Bulstoer(y,x1,x2,n,F) — дает более точное решение (методом Bulirsch-Stoer) Агрумкнты вышеуказанных функций:
y — вектор начальных условий
x1,x2 — границы интервала для поиска решения
n — количество точек на интервале
F(x,y) — вектор-функция первых производных При решении дифференциальных уравнений порядка выше первого (или систем уравнений, выше первого порядка) исходное уравнение (систему) необходимо преобразовать к системе дифференциальных уравнений первого порядка. В результате работы укзанных функций рассчитывается матрица, количество стобцов которой равно порядку уравнения +1(или сумме порядков уравнений в системе +1), а количество строк равно параметру n. Первый столбец содержит значения независимой переменной, второй — значение функции, третий — для диф. уравнений 2-го порядка — значение производной искомой функции (если решается система двух уравнений 1-го порядка, то третий столбец будет содержать значения второй функции). Для выделения решений (функций или их производных) можно воспользоваться стандартным оператором вывода столбцов матрицы M &lt &gt Если матрица правых частей дифференциальных уравнений почти вырождена, то такие системы называются жесткими. В этом случае решения, возвращаемые функцией rkfixed будет неустойчивым и для решения таких систем необходимо применять функции Stiffb , Stiffr Stiffb(y,x1,x2,n,F,J) — ищет решение диф. уравнения или системы дифференциальных уравнений методом Bulirsch-Stoer Stiffr(y,x1,x2,n,F,J) — ищет решение диф. уравнения или системы дифференциальных уравнений методом Rosenbrock

Первые пять аргументов такие же,как и при решении хорошо обусловленных систем дифференциальных уравнений . Дополнительный аргумент — матрица J размером nx(n+1), первый столбец которой содержит частные производные dF/dx, остальные столбцы и строки представляют собой матрицу Якоби dF/dy Пример решения жесткой системы дифференциальных уравнений. Для отыскания решения системы диф. уравнений только в конечной точке используются функции bulstoer,rkadapt, stiffb, stiffr (начинаются с прописной буквы).

Набор парамтров для этих функций :
bulstoer(y,x1,x2,acc,F,kmax,save)
rkadapt(y,x1,x2,acc,F,kmax,save)
stiffb(y,x1,x2,acc,F,J,kmax,save)
stiffr(y,x1,x2,acc,F,J,kmax,save) Первые три параметра и пятый (F) этих функций те же, что идля функции Rkadapt. Дополнительные параметры:
acc — параметр, контролирующий точность решения (реком. асс=0.001)
kmax — максимальное число промежуточных точек в которых ищется решение
save — минимально допустимый интервал между точками, в которых ищется решение

Решение граничных задач для обыкновенных дифференциальных уравнений.

Если для дифференциального уравнения n-го порядка k граничных условий заданы в начальной точке х1, а (n-k) граничных условий — в конечной точке х2, то такая задача называется краевой. В MathCAD реализованы две функции, позволяющие численно найти недостающие условия в точках х1 и х2. Двухточечная краевая задача Задача решается в два этапа. Сначала с помощью функции sbval находятся недостающие начальные значения, а затем применяется одна из выше описанных функций для решения стандартной задачи Коши на отрезке. sbval(v,x1,x2,F,load,score) — ищет недостающие начальные условия в точке х1
v — вектор началных приближений для искомых начальных значений в точке х1,
х1,х2 — граничные точки интервала
F(x,y) — вектор-столбец из n элементов, содержит правые части дифференциальных уравнений
load(x1,v) — вектор-столбец из n элементов, содержит начальные значения в точке х1; некоторые из значений- константы, другие неизвестны и будут найдены в процессе решения.
score(x2,y) — вектор-столбец размерности вектора v, содержащий разность между начальным условием в точке х2 и значеием искомого решения в этой точке. Краевая задача с условиями внутри интервала. На первом этапе используется функция balfit(V1,V2,x1,x2,xf,F,load1,load2,score) — ищет недостающие начальные условия в точках х1 и х2, сшивая решения, выходящие из этих точек, в точке xf
V1,V2 — вектора началных приближений для искомых начальных значений в точках х1 и х2
х1,х2 — граничные точки интервала
load1(x1,V1) — вектор-столбец из n элементов, содержит начальные значения в точке х1; некоторые из значений- константы, другие неизвестны и будут найдены в процессе решения
load2(x2,V2) — вектор-столбец из n элементов, содержит начальные значения в точке х2; некоторые из значений- константы, другие неизвестны и будут найдены в процессе решения.
score(xf,y) — вектор-столбец размерности n, содержащий разность между решениями, начинающимися в точках х1 и х2, в точке xf

9. РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В СРЕДЕ MATHCAD

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

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

Имеются два типа задач, которые возможно решать с помощью MathCAD:

— задачи Коши – для которых заданы значения функций в начальной точке интервала интегрирования уравнения;

— краевые задачи для которых заданы определенные соотношения сразу на обоих границах интервала.

MathCAD умеет решать только такие системы дифференциальных уравнений, которые могут быть представлены в стандартной форме (форме Коши). Для неизвестных функций y₀(t) ,y₁(t),…..y_N-1(t) система ОДУ должна быть записана в форме:

что эквивалентно следующему векторному представлению:

где Y и Y’ – соответствующие неизвестные векторные функции переменой t размера Nx1, F – векторная функция того же размера и количества переменных (N+1) (N компонент вектора и, возможно, t). Именно векторное представление (2) используется для ввода системы ОДУ в среде MathCAD.

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

где В – вектор начальных условий размера Nx1, составленный из y_i(t0).

Стандартные процедуры MathCAD применимы для систем ОДУ первого порядка, записанных в форме (2)–(3).

Дифференциальное уравнение N-го порядка

Для решения ОДУ порядка N≥1 в MathCAD предусмотрены две возможности:

— вычислительный блок Given/Odesolve – решение имеет вид функции от t;

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

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

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

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

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

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

На запись ОДУ в пределах вычислительного блока накладывается несколько

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

у(t0) = b или у ( N) (t0) = b.

Рис. 9.1. Пример решения задачи Коши

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

Пример решения задачи Коши для ОДУ второго порядка – системы генератора с затуханием (с параметром квадратичной нелинейности γ), приведен на рис. 9,1.

Символ производной допускается вводить как средствами панели Calculus (Вычисления), так и в виде штриха, набрав его с помощью сочетания клавиш +.

Результатом применения блока Given/Odesolve является функция y(t), определенная на промежутке (t0, t1). Следует воспользоваться обычными средствами MathCAD, чтобы построить ее график или получить значение функции в какой-либо точке указанного интервала.

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

Система N дифференциальных уравнений

При помощи MathCAD можно решать системы N>1 ОДУ первого порядка, если они записаны в стандартной форме (Коши) в виде векторного соотношении Y'(t)=F,(Y(t), t).

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

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

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

Bulstoer(y0, t0, t1, M, D) –метод Булирша–Штера:

у0 – вектор начальных значений в точке to размера Nx1;

t0 – начальная точка расчета;

t1 – конечная точка расчета;

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

D – векторная функция размера Nx1 двух аргументов – скалярного t и векторного у. При этом у – искомая векторная функция аргумента t того же размера Nx1.

Необходимо соблюдать регистр первой буквы рассматриваемых функций, поскольку это влияет на выбор алгоритма счета, в отличие от многих других встроенных функций MathCAD, например, Find ≡ find.

Каждая из приведенных функций выдает решение в виде матрицы размера (М+1)x(N+1). В ее левом столбце находятся значения аргумента t, делящие интервал на равномерные шаги, а в остальных N столбцах – значения искомых функций y₀(t), y₁(t), …, y_N-1(t), рассчитанные для этих значений аргумента. Поскольку всего точек

(помимо начальной) М, то строк в матрице решения будет всего М+1.

На рис. 56 приведен пример решения системы ОДУ генератора с затуханием при помощи функции rkfixed, результат расчетов представлен как содержимое матрицы и в виде графика. Точки, в которых получено решение, отмечены на графике кружками. Чтобы использовать другой численным алгоритм, достаточно поменять имя функции rkfixed в последней строке программы на другое (на практике более эффективны функции Rkadapt и Bulstoer).

Первая строка листинга представляет задание параметров модели, вторая – начального условия задачи Коши, а в третьей строке листинга определено число шагов, на которых рассчитывается решение. Самая важная – это предпоследняя строка листинга, в которой, собственно, определяется система ОДУ. Последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (0,40).

Функция D, входящая в число параметров встроенных функций для решения ОДУ, должна быть функцией обязательно двух аргументов. Второй ее аргумент должен быть вектором того же размера, что и сама функция D. Точно такой же размер должен быть и у вектора начальных значений у0. Векторную функцию D(t,y) следует определять через компоненты вектора у с помощью кнопки Subscript (Нижний индекс) с наборной панели Calculator (Калькулятор) или нажатием клавиши <[>.

Рис. 9.2. Пример решения

Размер полученной матрицы, представляющей решение, равен (M+1)x(N+1), т.е. 51х3. Просмотреть все компоненты матрицы u, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки.

Рис. 9.3. Уравнение

Решение одного уравнения (N=1)

Вычислительный блок Given/Odesolve выигрывает в простоте и в наглядности, однако иногда предпочтительнее решать ОДУ первого порядка с помощью встроенных функций rkfixed, Rkadapt или Bulstoer, например, при следующих обстоятельствах:

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

— ответ предпочтительнее получить в виде вектора, а не функции и т.п.

При решении одного ОДУ и само уравнение, и начальное условие можно задавать в скалярной форме (рис.

57). Результат выдается в виде матрицы размерности Мх2, которая состоит из двух столбцов: в одном находятся значения аргумента t (от t0 до t1 включительно), в другом – соответствующие значения искомой функции у(t).

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

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

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

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

bulstoer(у0, t0, t1, асс, D, k, s) – метод Булирша–Штера:

у0 – вектор начальных значений в точке to;

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

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

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

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

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

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

Рис. 9.4. Тестовый расчет

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

переменной М, в этой же строке оно выведено на экран.

В предпоследней программе осуществлен вывод решения системы ОДУ на конце интервала, т. е. в точке t=50 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующим местом вектора из предыдущей строки). Чтобы попробовать альтернативный численный метод, достаточно заменить имя функции bulstoer на rkadapt.

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

О численных методах

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

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

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

для решения единственного уравнения (любого порядка) используйте

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

— для стандартных нежестких систем используйте алгоритм Булирша– Штера (Bulstoer);

— для систем с участками быстро и медленно меняющихся решений используйте адаптивный алгоритм Рунге–Кутта (Rkadapt);

— в учебных целях и для решения несложных задач можно использовать алгоритм Рунге–Кутта с фиксированным шагом (Rkfixed);

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

Тем не менее, в MathCAD можно реализовать любой численный метод. Обычно решение заключается в нахождении ряда значений х_i и у_i искомой зависимости у(х) при i, изменяющемся от 0 до N при шаге изменения х, равном h. Рассмотрим три наиболее распространенных способа решения дифференциального уравнения, при которых h = const. Значения х₀ и у₀ должны быть известны как начальные условия, без них невозможно единственное решение.

Аналитическое решение обыкновенного дифференциального уравнения имеет вид: у(х)=А·ехр(х 2 /2). Погрешность метода определяется выражением: h n+1 , где n – порядок метода. Простой метод Эйлера имеет первый порядок, при этом погрешность метода составляет h 2 . Будучи самым простым из известных, этот метод имеет и наибольшую погрешность. При малых h и достаточно гладких решениях этот метод вполне приемлем для решения многих практических задач.

Простой метод Эйлера реализуется применением на каждом шаге вычислений следующих итерационных выражений:

Для уменьшения погрешности решения следует применять методы более высокого порядка. Модифицированный метод Эйлера является методом второго порядка, и его погрешность пропорциональна h 3 .

Данный метод реализуется следующими итерационными выражениями на каждом шаге вычислений:

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

Решение дифференциального уравнения второго порядка рассмотрим методом Рунге-Кутта четвертого порядка. Итерационные уравнения имеют вид:

Решить систему из двух дифференциальных уравнений простым методом Эйлера можно по итерационным уравнениям в векторной форме:

Модифицированный метод Эйлера имеет следующие уравнения:

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