Документация

Шаг 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 матрица. Матрица имеет 1 s и 0 s, где 1 s представляет вхождение переменной состояния. Например, 1 в положении, (2,3) означает второе уравнение, содержит третью переменную состояния T(t).

M = incidenceMatrix(eqns,vars)

M = 3×3

     1     0     1
     0     1     1
     1     1     0

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

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

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

Уменьшайте систему до системы первого порядка при помощи reduceDifferentialOrder. Функция 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, как описано в рабочем процессе, Решают Полулинейную Систему ДАУ.

Уменьшайте дифференциальный индекс ДАУ, описанных 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 ошибка или предупреждение, используйте альтернативный рабочий процесс, описанный в, Решают Полулинейную Систему ДАУ.

Часто, 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, смотрите, Решают ДАУ Используя Решатели Большой матрицы и Выбирают ODE Solver (MATLAB).

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: найдите начальные условия для решателей

Решатель ode15i требует начальных значений для всех переменных в указателе на функцию. Найдите начальные значения, которые удовлетворяют уравнения при помощи функции decic MATLAB. 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