Решение уравнений методом рунге кутта маткад

методами Эйлера и Рунге-Кутта в системе MathCAD

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

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

, .

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

,

.

Рассмотрим реализацию метода в MathCADна примере уравнения:

, , .

Аналитическое решение известно и имеет вид:

.

Краткие сведения о составлении программ в MathCAD

знак присваивает функции или переменной (они помещаются слева) выражение или число, которые помещаются справа. Набирается клавишей двоеточие «:» или из меню по цепочке View→ Toolbars→ Calculator.

знак обозначает последовательное изменение переменной через единицу от значения слева до значения справа. Набирается клавишей точка с запятой «;».

матрица вставляется командой меню Insert→ Matrixили клавишами Ctrl-M. Нижний индекс добавляется клавишей квадратная скобка «[».

Графиквставляется командой Insert→ Graph→ X-YPlotили клавишей «@».

Для удобства в работе рекомендуется отключить автоматическое вычисление, убрав галочку с опции меню Tools→ Calculate→ AutomaticCalculation. Тогда расчет не будет выполняться в ходе набора программы, а запуститься только после нажатия кнопки Calculate, расположенной на панели инструментов (в виде значка ).

Ниже приведена расчетная программа. Повторите её. Получите графики с тем же форматом линий. Формат линий графика можно изменить, открыв с помощью ПК мыши контекстное меню и выбрав Format… → Traces.

Программа для простого метода Эйлера

Следующая программа реализует модифицированный метод Эйлера. Отличие от простого метода заключается в итерационных уравнениях.

Программа для модифицированного метода Эйлера

Метод Рунге-Кутта 4-го порядка используется в тех случаях, когда необходима высокая точность расчетов, недостигаемая методами Эйлера.

Программа для метода Рунге-Кутта

Решение дифференциальных уравнений 2-го порядкаметодом Рунге-Кутта.

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

, , .

Преобразуем задачу к системе из двух уравнений:

, ,

, .

Тогда получим следующее обобщение итерационной схемы:

. ,

, ,

, ,

, ,

, .

Отметим, что значения на каждом следующем шаге рассчитываются по значениям, полученным на предыдущем. Кроме того, использованы прежние правила «взвешивания» коэффициентов при усреднении.

Пример математической модели с дифференциальным уравнением 2-го порядка

Рассмотрим уравнение колебательного процесса при наличии внешнего периодического воздействия:

,

где t– время, и искомой является зависимость ;

– круговая частота собственных колебаний;

– круговая частота внешнего воздействия с амплитудой «a».

Если , то общее решение уравнения имеет вид (проверьте подстановкой):

,

где Aи – произвольные постоянные. Частное решение выбирается заданием значений этих постоянных. Второе слагаемое решения показывает, что с течением времени амплитуда колебаний неограниченно возрастает. Это явление называется резонансом.

Когда , общее решение имеет вид:

.

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

Моделирование резонансных колебаний

Методом Рунге-Кутта найдем решение задачи:

, , , .

Согласно изложенной выше теории, аналитическое решение уравнения имеет вид:

.

Ниже приведен алгоритм расчета и его реализация в MathCAD.

Программа расчета резонансных колебаний методом Рунге-Кутта

Задание для самостоятельного выполнения

Найти решение уравнения вынужденных колебаний:

, , , .

Решение представить в виде графика. Для сравнения привести и график точного решения (также как это было сделано для резонансных колебаний).

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

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

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

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

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

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

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

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

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

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

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

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

При решении дифференциального уравнения первого порядка нужно создать вектор начальных условий из одного элемента Y1, который затем используется при формировании вектора-функции правой части дифференциального уравнения. При обращении к функции 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. Решение дифференциальных уравнений и систем в 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.

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