Шаг 1: Задайте уравнения и переменные

Следующий рисунок показывает рабочий процесс ДАУ путем решения ДАУ для маятника.

Переменные состояния:

Переменные:

Система уравнений ДАУ:

Задайте независимые переменные и переменные состояния при помощи 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);

Шаг 2: Уменьшите дифференциальный порядок

2.1 (необязательно) Проверяйте частоту переменных

Этот шаг является необязательным. Можно проверить, где переменные происходят в системе ДАУ, просмотрев матрицу заболеваемости. Этот шаг находит любые переменные, которые не происходят в вашем входе и могут быть удалены из vars вектор.

Отобразите матрицу инцидентности при помощи incidenceMatrix. Область выхода incidenceMatrix имеет строку для каждого уравнения и столбец для каждой переменной. Поскольку система имеет три уравнения и три переменные состояния, incidenceMatrix возвращает 3-by- 3 матрица. Матрица имеет 1s и 0s, где 1s представляет вхождение переменной состояния. Для примера, 1 в позиционном (2,3) означает, что второе уравнение содержит третью переменную состояния T(t).

M = incidenceMatrix(eqns,vars)

M = 3×3

     1     0     1
     0     1     1
     1     1     0

Если столбец матрицы инцидентности все 0s, тогда эта переменная состояния не происходит в системе ДАУ и должна быть удалена.

2.2 Уменьшите дифференциальный порядок

Дифференциальный порядок системы ДАУ является самым высоким дифференциальным порядком ее уравнений. Чтобы решить ДАУ с помощью MATLAB, дифференциальный порядок должен быть уменьшен до 1. Здесь первое и второе уравнения имеют производные второго порядка x(t) и y(t). Таким образом, дифференциальный порядок 2.

Сведите систему к системе первого порядка при помощи reduceDifferentialOrder. The reduceDifferentialOrder функция заменяет производные новыми переменными, такими как Dxt(t) и Dyt(t). Правая сторона выражений в eqns является 0.

[eqns,vars] = reduceDifferentialOrder(eqns,vars)

eqns = 
$$\left(\begin{array}{c}m\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dxt}\left(t\right)-\frac{T\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ g\hspace{0.17em}m+m\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)-\frac{T\left(t\right)\hspace{0.17em}y\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ \mathrm{Dxt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }x\left(t\right)\\ \mathrm{Dyt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

vars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ T\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\end{array}\right)$$

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

3.1 Проверяйте дифференциальный индекс системы

Проверяйте дифференциальный индекс системы ДАУ при помощи isLowIndexDAE. Если индекс 0 или 1, затем isLowIndexDAE возвращает логический 1 (true) и можно пропустить шаг 3.2 и перейти к шагу 4. Преобразуйте системы ДАУ в указатели на функции MATLAB. Здесь, isLowIndexDAE возвращает логический 0 (false), что означает, что дифференциальный индекс больше 1 и должна быть уменьшена.

isLowIndexDAE(eqns,vars)

ans = logical
   0

3.2 Уменьшение дифференциального индекса с помощью reduceDAEIndex

Чтобы уменьшить дифференциальный индекс, reduceDAEIndex функция добавляет новые уравнения, которые получают из входных уравнений, а затем заменяет производные более высокого порядка новыми переменными. Если reduceDAEIndex отказывает и выдает предупреждение, затем использует альтернативную функцию reduceDAEToODE как описано в рабочем процессе Solve Semilinear DAE System.

Уменьшите дифференциальный индекс ДАУ, описываемых eqns и vars.

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)

DAEs = 
$$\left(\begin{array}{c}m\hspace{0.17em}\mathrm{Dxtt}\left(t\right)-\frac{T\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ g\hspace{0.17em}m+m\hspace{0.17em}\mathrm{Dytt}\left(t\right)-\frac{T\left(t\right)\hspace{0.17em}y\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ \mathrm{Dxt}\left(t\right)-{\mathrm{Dxt}}_1\left(t\right)\\ \mathrm{Dyt}\left(t\right)-{\mathrm{Dyt}}_1\left(t\right)\\ 2\hspace{0.17em}{\mathrm{Dxt}}_1\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}{\mathrm{Dyt}}_1\left(t\right)\hspace{0.17em}y\left(t\right)\\ 2\hspace{0.17em}y\left(t\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }{\mathrm{Dyt}}_1\left(t\right)+2\hspace{0.17em}{{\mathrm{Dxt}}_1\left(t\right)}^2+2\hspace{0.17em}{{\mathrm{Dyt}}_1\left(t\right)}^2+2\hspace{0.17em}\mathrm{Dxt1t}\left(t\right)\hspace{0.17em}x\left(t\right)\\ \mathrm{Dxtt}\left(t\right)-\mathrm{Dxt1t}\left(t\right)\\ \mathrm{Dytt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }{\mathrm{Dyt}}_1\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

DAEvars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ T\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\\ \mathrm{Dytt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)\\ {\mathrm{Dxt}}_1\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)\\ \mathrm{Dxt1t}\left(t\right)\end{array}\right)$$

Если reduceDAEIndex выдает ошибку или предупреждение, используйте альтернативный рабочий процесс, описанный в Solve Semilinear DAE System.

Часто, reduceDAEIndex представляет избыточные уравнения и переменные, которые можно исключить. Исключить избыточные уравнения и переменные можно используя команду reduceRedundancies.

[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)

DAEs = 
$$\left(\begin{array}{c}-\frac{T\left(t\right)\hspace{0.17em}x\left(t\right)-m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxtt}\left(t\right)}{r}\\ \frac{g\hspace{0.17em}m\hspace{0.17em}r-T\left(t\right)\hspace{0.17em}y\left(t\right)+m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dytt}\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ 2\hspace{0.17em}\mathrm{Dxt}\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}\mathrm{Dyt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ 2\hspace{0.17em}{\mathrm{Dxt}\left(t\right)}^2+2\hspace{0.17em}{\mathrm{Dyt}\left(t\right)}^2+2\hspace{0.17em}\mathrm{Dxtt}\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}\mathrm{Dytt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ \mathrm{Dytt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

DAEvars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ T\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\\ \mathrm{Dytt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)\end{array}\right)$$

Проверьте дифференциальный индекс новой системы. Теперь, isLowIndexDAE возвращает логический 1 (true), что означает, что дифференциальный индекс системы 0 или 1.

isLowIndexDAE(DAEs,DAEvars)

ans = logical
   1

Шаг 4: Преобразование систем ДАУ в указатели на функции MATLAB

Этот шаг создает указатели на функцию для решателя MATLAB ® ODE ode15i, который является решателем общего назначения. Чтобы использовать специализированные большие матрицы, такие как ode15s и ode23t, см. Решение ДАУ Использование решателей Большой матрицы и Выбор решателя ОДУ.

reduceDAEIndex выводит вектор уравнений в DAEs и вектор переменных в DAEvars. Как использовать ode15i, вам нужен указатель на функцию, который описывает систему ДАУ.

Во-первых, уравнения в DAEs может содержать символьные параметры, которые не заданы в векторе переменных DAEvars. Найдите эти параметры при помощи setdiff на выходе symvar от DAEs и DAEvars.

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars)

extraParams = $$\left(\begin{array}{ccc}g& m& r\end{array}\right)$$

Дополнительные параметры, которые вы должны задать, это масса m, радиус r, и ускорение свободного падения g.

Создайте указатель на функцию при помощи daeFunction. Задайте дополнительные символьные параметры как дополнительные входные параметры daeFunction.

f = daeFunction(DAEs,DAEvars,g,m,r);

Остальная часть рабочего процесса является чисто числовой. Установите значения параметров и создайте указатель на функцию для ode15i.

g = 9.81;
m = 1;
r = 1;
F = @(t,Y,YP) f(t,Y,YP,g,m,r);

Шаг 5: Найти начальные условия для решателей

The ode15i решатель требует начальных значений для всех переменных в указателе на функцию. Найдите начальные значения, которые удовлетворяют уравнениям при помощи MATLAB decic функция. decic принимает догадки (которые могут не удовлетворять уравнениям) для начальных условий и пытается найти удовлетворительные начальные условия, используя эти догадки. decic может потерпеть неудачу, в этом случае необходимо вручную задать согласованные начальные значения для вашей задачи.

Сначала проверьте переменные в DAEvars.

DAEvars

DAEvars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ T\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\\ \mathrm{Dytt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)\end{array}\right)$$

Здесь, Dxt(t) является первой производной x(t), Dytt(t) является второй производной y(t)и так далее. В 7 7 переменных-by- 1 вектор. Поэтому догадки для начальных значений переменных и их производных также должны быть 7-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; 0; 0];
yp0est = zeros(7,1);

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

opt = odeset('RelTol', 10.0^(-7),'AbsTol',10.0^(-7));

Найдите непротиворечивые начальные значения для переменных и их производных при помощи decic.

[y0,yp0] = decic(F,0,y0est,[],yp0est,[],opt)

y0 = 7×1

    0.4771
   -0.8788
   -8.6214
         0
    0.0000
   -2.2333
   -4.1135

yp0 = 7×1

         0
    0.0000
         0
         0
   -2.2333
         0
         0

Шаг 6: Решить ДАУ используя ode15i

Решите систему, интегрирующуюся за временной промежуток 0 ≤ t ≤ 0.5. Добавить на график линии сетки и легенду.

[tSol,ySol] = ode15i(F,[0 0.5],y0,yp0,opt);
plot(tSol,ySol(:,1:origVars),'LineWidth',2)

for k = 1:origVars
  S{k} = char(DAEvars(k));
end

legend(S,'Location','Best')
grid on

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

Задайте r на 2 и регенерируйте указатель на функцию и начальные условия.

r = 2;
F = @(t,Y,YP)f(t,Y,YP,g,m,r);

y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0; 0; 0];
[y0,yp0] = decic(F,0,y0est,[],yp0est,[],opt);

Решить систему для нового значения параметров.

[tSol,y] = ode15i(F,[0 0.5],y0,yp0,opt);
plot(tSol,y(:,1:origVars),'LineWidth',2)

for k = 1:origVars
  S{k} = char(DAEvars(k));
end
legend(S,'Location','Best')
grid on