Шаг 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