Функции mathcad для решения обыкновенных дифференциальных уравнений

Функции mathcad для решения обыкновенных дифференциальных уравнений

Электронный курс по 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

Функции mathcad для решения обыкновенных дифференциальных уравнений

Для решения дифференциальных уравнений Mathcad предоставляет пользователю библиотеку встроенных функций Differential Equation Solving, предназначенных для численного решения дифференциальных уравнений.

ODESOLVE

Встроенная функция odesolve предназначена для решения дифференциальных уравнений, линейных относительно старшей производной. В отличие от других функций библиотеки Differential Equation Solving, odesolve решает дифференциальные уравнения, записанные в общепринятом в математической литературе виде.

• Функция odesolve решает для уравнений вида
  a(x) y(n) + F(x, y, y’ , . y (n-1) )=f(x)
  задачу Коши
  y(x₀ )=y₀ , y'(x₀ )=y_0,1 , y»(x₀ )=y_0,2 , . y (n-1) (x₀ )=y_0,n-1
  или простейшую граничную задачу
  y (k) (a)=y_a,k , y (m) (b)=y_b,k , 0 + ), а для записи производных можно использовать как оператор дифференцирования, так и знак производной, например, вторую производную можно вводить в виде или в виде y»(x). При этом необходимо обязательно записывать аргумент искомой функции.
• Для того чтобы вывести в рабочий документ значения решения в любой точке промежутка интегрирования, достаточно ввести имя функции Y, указать в скобках значение аргумента и знак равенства.
• Значения решения в любой точке промежутка интегрирования можно использовать в дальнейших вычислениях, достаточно ввести в нужном месте имя функции Y, указав в скобках значение аргумента.

Полную информацию о правилах использования функции odesolve можно получить во встроенном справочнике Mathcad в разделе Overview fnd Tutorials.

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

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

Рассмотрим задачу Коши:

Численное решение этой задачи состоит в построении таблицы приближенных значений
y_i,1 , y_i,2 , . y_i,N
решения y₁ (x), y₂ (x), . y_N (x)
на отрезке [x₀ , x_N ] в точках
x₁ , x₂ , . x_N, которые называются узлами сетки.
Обозначив

, ,

,

,

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

, .

В Mathcad решить задачу Коши для такой системы можно с помощью следующих функций:

• rkfixed(y, x1, x2, npoints, D) —решение задачи на отрезке методом Рунге—Кутты с постоянным шагом;
• Rkadapt(y, x1, x2, npoints, D) —решение задачи на отрезке методом Рунге—Кутты с автоматическим выбором шага;
• rkadapt(y, x1, x2, acc, npoints, D, kmax, save) —решения задачи в заданной точке методом Рунге-Кутты с автоматическим выбором шага;
• Bulstoer(y, x1, x2, npoints, D) —решение задачи на отрезке методом Булирша-Штера;
• bulstoer(y, x1, x2, acc, npoints, D, kmax, save) —решение задачи в заданной точке методом Булирша—Штера;
• Stiffr(y, x1, x2, acc, D, J) — решение задачи для жестких систем на отрезке с использованием алгоритма Розенброка;
• stiffr(y, x1, x2, acc, D, J, kmax, save) —решения задач для жестких систем на отрезке с использованием алгоритма Розенброка;
• Stiffb(y, x1, x2, acc, D, J) —решение задачи для жестких систем на отрезке с использованием алгоритма Булирша—Штера;
• stiffb(y, x1, x2, acc, D, J, kmax, save) —решение задач для жестких систем в заданной точке с использованием алгоритма Булирша—Штера.

Смысл параметров для всех функций одинаков и определяется математической постановкой задачи:
y — вектор начальных условий , ;
x1, x2 — начальная и конечная точки отрезка интегрирования системы; для функций, вычисляющих решение в заданной точке, x1 — начальная точка, x2 — заданная точка;
npoints — число узлов на отрезке [x1, x]; при решении задачи на отрезке результат содержит npoints+1 строку;
D — имя вектор-функции D(x,y) правых частей , ; ( имя D – от Derivative — производная, имя вектора, содержащего выражения для производных (derivatives) искомого решения);
J — имя матрицы-функции J(x,y) размерности n x (n+1), в первом столбце которой хранятся выражения частных производных по x правых частей системы, а в остальных n столбцах содержится матрица Якоби правых частей:
.
acc — параметр, контролирующий погрешность решения при автоматическом выборе шага интегрирования (если погрешность решения больше acc, то шаг сетки уменьшается; шаг уменьшается до тех пор, пока его значение не станет меньше save );
kmax — максимальное число узлов сетки, в которых может быть вычислено решение задачи на отрезке, максимальное число строк в результате;
save — наименьшее допустимое значение шага неравномерной сетки.

Результат работы функции — матрица, содержащая n+1; ее первый столбец содержит координаты узлов сетки, второй столбец — вычисленные приближенные значения решения y₁ (x) в узлах сетки, (k+1) -й — значения решенияy_k (x) в узлах сетки.

При решении задачи Коши для дифференциального уравнения первого порядка результат вычислений всех приведенных выше функций — матрица, в первом столбце которой содержатся координаты узлов сетки x₀ , x₁ , . x_N, а во втором — значения приближенного решения в соответствующих узлах.

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

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

Запишем автономную систему второго порядка

.

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

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

Тема 7. Решение дифференциальных уравнений и систем в MathCad

Краткие теоретические сведения

Для решения дифференциальных уравнений с начальными условиями система Mathcad имеет ряд встроенных функций:

rkfixed – функция для решения ОДУ и систем ОДУ методом Рунге–Кутта четвертого порядка с постоянным шагом;

Rkadapt – функция решения ОДУ и систем ОДУ методом Рунге–Кутта с переменным шагом;

Odesolve – функция, решающая ОДУ блочным методом.

Ниже приведено описание стандартной функции rkfixed с указанием параметров функции.

y – вектор начальных условий из k элементов ( k – количество уравнений в системе);

x1 и x2 – левая и правая границы интервала, на котором ищется решение ОДУ или системы ОДУ;

p – число точек внутри интервала (x1, x2), в которых ищется решение;

D – вектор, состоящий из k-элементов, который содержит первую производную искомой функции или первые производные искомых функций, если речь идет о решении системы.

Результатом работы функции является матрица из p +1 строк, первый столбец которой содержит точки, в которых получено решение, а остальные столбцы – сами решения.

На рисунке 2.7.1 приведены конкретные примеры решения различных дифференциальных уравнений и систем ОДУ в MathCAD .

Рисунок 2.7.1 – Примеры решения дифференциальных уравнений и систем

При решении дифференциального уравнения первого порядка нужно создать вектор начальных условий из одного элемента Y 1 , который затем используется при формировании вектора-функции правой части дифференциального уравнения. При обращении к функции rkfixed указывается имя вектора Y , границы интервала, на котором ищется решение уравнения, например, (0 ; 2), количество точек, в которых ищется решение – 100, вектор-функция, описывающая правую часть дифференциального уравнения – D . В результате получается матрица z , в первом столбце которой содержатся значения аргумента искомой функции, во втором – значения самой результирующей функции. При построении графика функции первый столбец полученной матрицы указывается как аргумент, второй столбец – как функция.

При решении системы дифференциальных уравнений нужно создать вектор начальных условий из двух элементов, например, вектор v , который затем используется при формировании вектора-функции правой части дифференциального уравнения. При обращении к функции rkfixed указывается имя вектора v , и границы интервала, на котором ищется решение уравнения, например, (0 ; 5), количество точек, в которых ищется решение – 100, вектор-функция, описывающая правую часть дифференциального уравнения – D . В результате получается матрица s , в первом столбце которой содержатся значения аргумента искомых функций, во втором и третьем столбцах – значения самих функций при соответствующем значении аргумента. При построении графика можно воспользоваться первым столбцом полученной матрицы как аргументом, а вторым и третьим столбцами – как функциями.

На рисунке 2.7.2 приведен пример решения дифференциального уравнения второго порядка с использованием функции rkfixed . Необходимо решить дифференциальное уравнение второго порядка с заданными начальными условиями вида:

Рисунок 2.7.2 – Пример решения дифференциальных уравнений второго порядка с помощью rkfixed

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

Документ формируется точно так же, как и при решении системы ОДУ.

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

Практическая часть темы 7

7.1 Решение дифференциальных уравнений первого порядка

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

q сформировать вектор начальных условий из одного элемента, присвоив начальное значение искомой функции переменной с индексом, например: или (в зависимости от значения переменной ORIGIN );

q определить вектор-функцию из одного элемента, которая содержит первую производную неизвестной функции:

· набрать имя функции с двумя параметрами: первый параметр – аргумент искомой функции (независимая переменная), второй – имя вектора, содержащего искомую функцию (можно использовать имя вектора начальных условий), например, D ( x , Y );

· набрать оператор «:=» и выражение для первой производной (выразить из дифференциального уравнения), в котором вместо имени искомой функции подставлен первый элемент вектора-параметра, например, для уравнения вектор-функция будет определятся следующим образом: ( если ORIGIN = 0 , подставлять );

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

· первый – имя вектора начальных условий,

· второй – левая граница интервала, на котором ищется решение, в виде числовой константы,

· третий – правая граница интервала, на котором ищется решение, в виде числовой константы,

· четвертый – количество точек, в которых ищется решение,

· пятый – имя вектора-функции, описывающего первую производную, без параметров;

например: ,

(в результате получится матрица Z , в первом столбце которой содержатся значения аргумента искомой функции, во втором – значения самой функции);

q вывести матрицу, содержащую решение ДУ с помощь оператора «=», например: Z = ;

q построить график найденной функции ( см. тему 5 ), указав в качестве аргумента по оси абсцисс столбец , а в качестве значения функции по оси ординат – столбец ( если ORIGIN = 0 , набирать соответственно и ).

Пример 7.1 Найти численное решение дифференциального уравнения первого порядка на интервале от 0.2 до 5 в 1000 точках, при начальном условии y (0)=0.1.

Выполнить графическую интерпретацию результатов.

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

Последовательность действий для р ешения системы дифференциальных уравнений первого порядка такова (описана для значения ORIGIN =0 ):

q перейти в исходной системе уравнений к однотипным обозначениям функций и выразить первые производные,

например, систему можно преобразовать в ;

q в документе MathCad сформировать вектор начальных условий, количество элементов которого равно количеству уравнений системы, присвоив его некоторой переменной (см. тему 2);

например, ;

q определить вектор-функцию, которая содержит первые производные искомых функций:

· набрать имя функции с двумя параметрами: первый параметр – аргумент искомых функций (независимая переменная), второй – имя вектора, содержащего искомые функции (можно использовать имя вектора начальных условий), например, D ( t , V );

(Замечание: если независимая переменная явно не присутствует в системе, то в качестве ее имени можно выбрать любую переменную)

· набрать оператор «:=» и вставить шаблон вектора, количество элементов которого равно количеству уравнений системы (см. тему 2)

· набрать в качестве элементов вектора правые части системы уравнений, в которых искомые функции представлены соответствующими элементами вектора-параметра, например,

;

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

· первый – имя вектора начальных условий,

· второй – левая граница интервала, на котором ищется решение, в виде числовой константы,

· третий – правая граница интервала, на котором ищется решение, в виде числовой константы,

· четвертый – количество точек, в которых ищется решение,

· пятый – имя вектора-функции, описывающего первые производные, без параметров;

например: ,

(в результате получится матрица Z , в первом столбце которой содержатся значения аргумента искомых функций, во втором – значения первой функции, в третьем – значения второй функции и т. д.);

q вывести матрицу, содержащую решение системы ДУ с помощь оператора «=», например: Z = ;

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

Пример 7.2 Найти решение системы дифференциальных уравнений

на интервале от 0 до 0.5 в 1000 точках, при следующих начальных условиях: x (0)=0.1 и y (0)=1.

Выполнить графическую интерпретацию результатов.