Шаг 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- 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. 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 ложь), что означает, что дифференциальный индекс больше 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.

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 переменных- 1 вектор. Поэтому предположениями для начальных значений переменных и их производных должен также быть 7- 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