Этот рабочий процесс является альтернативным рабочим процессом для решения полулинейных дифференциальных алгебраических уравнений (ДАУ), используемых только если reduceDAEIndex
сбой в стандартном рабочем процессе с предупреждающим сообщением: The index of the reduced DAEs is larger than 1. [daetools::reduceDAEIndex]
. Для стандартного рабочего процесса см. Решение дифференциальных алгебраических уравнений (ДАУ).
Завершите шаги 1 и 2 в Solve Differential Algebraic Equations (ДАУ), прежде чем начинать другие шаги. Затем, на шаге 3, если reduceDAEIndex
fails, уменьшить дифференциальный индекс используя reduceDAEToODE
. Преимущество reduceDAEToODE
заключается в том, что она надежно снижает полунедельные ДАУ до ОДУ (ДАУ индекса 0
). Однако эта функция медленнее и работает только на полумесячных системах ДАУ. reduceDAEToODE
может завершиться неуспешно, если система не является полулинейной.
Чтобы решить систему ДАУ, выполните эти шаги.
Система уравнений ДАУ:
Задайте независимые переменные и переменные состояния при помощи syms
.
syms x(t) y(t) T(t) m r g
Задайте уравнения при помощи оператора = =.
eqn1 = m*diff(x(t), 2) == T(t)/r*x(t); eqn2 = m*diff(y(t), 2) == T(t)/r*y(t) - m*g; eqn3 = x(t)^2 + y(t)^2 == r^2; eqns = [eqn1 eqn2 eqn3];
Поместите переменные состояния в вектор-столбец. Сохраните количество исходных переменных для ссылки.
vars = [x(t); y(t); T(t)]; origVars = length(vars);
Дифференциальный порядок системы ДАУ является самым высоким дифференциальным порядком ее уравнений. Чтобы решить ДАУ с помощью MATLAB, дифференциальный порядок должен быть уменьшен до 1
. Здесь первое и второе уравнения имеют производные второго порядка x(t)
и y(t)
. Таким образом, дифференциальный порядок 2
.
Сведите систему к системе первого порядка при помощи reduceDifferentialOrder
. The reduceDifferentialOrder
функция заменяет производные новыми переменными, такими как Dxt(t)
и Dyt(t)
. Правая сторона выражений в eqns
является 0
.
[eqns,vars] = reduceDifferentialOrder(eqns,vars)
eqns =
vars =
reduceDAEToODE
Чтобы уменьшить дифференциальный индекс ДАУ, описываемых eqns
и vars
, использовать reduceDAEToODE
. Чтобы уменьшить индекс, reduceDAEToODE
добавляет новые переменные и уравнения в систему. reduceDAEToODE
также возвращает ограничения, которые являются условиями, которые помогают найти начальные значения, чтобы гарантировать, что получившиеся ОДУ равны исходным ДАУ.
[ODEs,constraints] = reduceDAEToODE(eqns,vars)
ODEs =
constraints =
Из выхода reduceDAEToODE
, у вас есть вектор уравнений в ODEs
и вектор переменных в vars
. Как использовать ode15s
или ode23t
, вам нужно два указателя на функцию: один, представляющая большая матрица ОДУ системы, и другой, представляющий вектор, содержащий правые стороны большой матрицы уравнений. Эти указатели на функцию являются эквивалентным матричным представлением системы ODE, где M (t, y (t)) y "(t) = f (t, y (t)).
Найдите эти указатели на функцию при помощи massMatrixForm
для получения большой матрицы massM
(M в уравнении) и правую сторону f
.
[massM,f] = massMatrixForm(ODEs,vars)
massM =
f =
Уравнения в ODEs
может содержать символьные параметры, которые не заданы в векторе переменных vars
. Найдите эти параметры при помощи setdiff
на выходе symvar
от ODEs
и vars
.
pODEs = symvar(ODEs); pvars = symvar(vars); extraParams = setdiff(pODEs, pvars)
extraParams =
Дополнительные параметры, которые вы должны задать, это масса m
, радиус r
, и ускорение свободного падения g
.
Преобразование massM
и f
для указателей на функцию с помощью odeFunction
. Задайте дополнительные символьные параметры в качестве дополнительных входов для odeFunction
.
massM = odeFunction(massM, vars, m, r, g); f = odeFunction(f, vars, m, r, g);
Остальная часть рабочего процесса является чисто числовой. Установите значения параметров и замените значения параметров в DAEs
и constraints
.
m = 1; r = 1; g = 9.81; ODEsNumeric = subs(ODEs); constraintsNumeric = subs(constraints);
Создайте указатель на функцию, подходящий для входа в ode15s
или ode23s
.
M = @(t,Y) massM(t,Y,m,r,g); F = @(t,Y) f(t,Y,m,r,g);
ode15s
и ode23t
Решатели требуют начальных значений для всех переменных в указателе на функцию. Найдите начальные значения, которые удовлетворяют уравнениям при помощи MATLAB ® decic
функция. The decic
принимает догадки (которые могут не удовлетворять уравнениям) для начальных условий и пытается найти удовлетворительные начальные условия, используя эти догадки. decic
может потерпеть неудачу, в этом случае необходимо вручную задать согласованные начальные значения для вашей задачи.
Сначала проверьте переменные в vars
.
vars
vars =
Здесь, Dxt(t)
является первой производной x(t)
и так далее. В 5
5 переменных-by-
1
вектор. Поэтому догадки для начальных значений переменных и их производных также должны быть 5
-by- 1
векторы.
Предположим, что начальное угловое перемещение маятника составляет 30 ° или pi/6
, и источник координат находится в точке подвеса маятника. Учитывая, что мы использовали радиус r
от 1
, начальное горизонтальное положение x(t)
является r*sin(pi/6)
. Начальное вертикальное положение y(t)
является -r*cos(pi/6)
. Задайте эти начальные значения переменных в векторе y0est
.
Произвольно установите начальные значения остальных переменных и их производных в 0
. Это плохие догадки. Однако их достаточно для этой проблемы. В вашей задаче, если decic
ошибки, затем предоставьте лучшие догадки и обратитесь к decic
страница.
y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0]; yp0est = zeros(5,1);
Создайте набор опций, содержащий большую матрицу M
системы и определяет числовые допуски для численного поиска.
opt = odeset('Mass', M, 'RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
Найти начальные значения, сопоставимые с системой ОДУ и с алгебраическими ограничениями, используя decic
. Значение параметра [1,0,0,0,1]
в этом вызове функции фиксирует первый и последний элемент в y0est
, так что decic
не изменяет их во время численного поиска. Здесь это исправление необходимо для обеспечения decic
находит удовлетворительные начальные условия.
[y0, yp0] = decic(ODEsNumeric, vars, constraintsNumeric, 0,...
y0est, [1,0,0,0,1], yp0est, opt)
y0 = 5×1
0.5000
-0.8660
-8.4957
0
0
yp0 = 5×1
0
0
0
-4.2479
-2.4525
Теперь создайте набор опций, содержащий большую матрицу M
системы и вектор yp0
допустимых начальных значений для производных. Вы будете использовать этот набор опций при решении системы.
opt = odeset(opt, 'InitialSlope', yp0);
ode15s
или ode23t
Решите систему, интегрирующуюся за временной промежуток 0
≤ t
≤ 0.5
. Добавить на график линии сетки и легенду. Использование ode23s
путем замены ode15s
с ode23s
.
[tSol,ySol] = ode15s(F, [0, 0.5], y0, opt); plot(tSol,ySol(:,1:origVars),'-o') for k = 1:origVars S{k} = char(vars(k)); end legend(S, 'Location', 'Best') grid on
Решить систему для различных значений параметров путем установки нового значения и регенерации указателя на функцию и начальных условий.
Задайте r
на 2
и повторите шаги.
r = 2; ODEsNumeric = subs(ODEs); constraintsNumeric = subs(constraints); M = @(t,Y) massM(t,Y,m,r,g); F = @(t,Y) f(t,Y,m,r,g); y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0]; opt = odeset('Mass', M, 'RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7)); [y0, yp0] = decic(ODEsNumeric, vars, constraintsNumeric, 0,... y0est, [1,0,0,0,1], yp0est, opt); opt = odeset(opt, 'InitialSlope', yp0);
Решить систему для нового значения параметров.
[tSol,ySol] = ode15s(F, [0, 0.5], y0, opt); plot(tSol,ySol(:,1:origVars),'-o') for k = 1:origVars S{k} = char(vars(k)); end legend(S, 'Location', 'Best') grid on
daeFunction
| decic
| findDecoupledBlocks
| incidenceMatrix
| isLowIndexDAE
| massMatrixForm
| odeFunction
| reduceDAEIndex
| reduceDAEToODE
| reduceDifferentialOrder
| reduceRedundancies