Решите полулинейную систему ДАУ

Этот рабочий процесс является альтернативным рабочим процессом к решению полулинейных ДАУ, используемых, только если reduceDAEIndex перестал работать в стандартном рабочем процессе с предупреждением ниже. Для стандартного рабочего процесса смотрите, Решают Дифференциальные Алгебраические уравнения (ДАУ).

Warning: The index of the reduced DAEs is larger...
than 1. [daetools::reduceDAEIndex]

Чтобы решить вашу систему ДАУ завершают эти шаги.

Шаг 1. Уменьшайте дифференциальный индекс с `reduceDAEToODE`

Полные шаги 1 и 2 в Решают Дифференциальные Алгебраические уравнения (ДАУ) прежде, чем начать этот шаг. Затем на шаге 3, когда reduceDAEIndex перестанет работать, уменьшайте дифференциальный индекс с помощью reduceDAEToODE. Преимущество reduceDAEToODE состоит в том, что он надежно уменьшает полулинейные ДАУ до ОДУ (ДАУ индекса 0). Однако эта функция медленнее и работает только над полулинейными системами ДАУ. reduceDAEToODE может перестать работать, если система не полулинейна.

Чтобы уменьшать дифференциальный индекс ДАУ, описанных eqns и vars, используйте reduceDAEToODE. Чтобы уменьшать индекс, reduceDAEToODE добавляет новые переменные и уравнения к системе. reduceDAEToODE также возвращает ограничения, которые являются условиями, что справка находит, что начальные значения гарантируют, что получившиеся ОДУ равны начальным ДАУ.

[ODEs,constraints] = reduceDAEToODE(eqns,vars)

ODEs =
                             Dxt(t) - diff(x(t), t)
                             Dyt(t) - diff(y(t), t)
                  m*diff(Dxt(t), t) - (T(t)*x(t))/r
          m*diff(Dyt(t), t) - (T(t)*y(t) - g*m*r)/r
        -(4*T(t)*y(t) - 2*g*m*r)*diff(y(t), t) -...
           diff(T(t), t)*(2*x(t)^2 + 2*y(t)^2) -...
                     4*T(t)*x(t)*diff(x(t), t) -...
                  4*m*r*Dxt(t)*diff(Dxt(t), t) -...
                       4*m*r*Dyt(t)*diff(Dyt(t), t)
 
constraints =
 2*g*m*r*y(t) - 2*T(t)*y(t)^2 - 2*m*r*Dxt(t)^2 -...
                     2*m*r*Dyt(t)^2 - 2*T(t)*x(t)^2
                              r^2 - y(t)^2 - x(t)^2
                       2*Dxt(t)*x(t) + 2*Dyt(t)*y(t)

Шаг 2. ОДУ к Указателям на функции для ode15s и ode23t

От вывода reduceDAEToODE у вас есть вектор уравнений в ODEs и вектор переменных в vars. Чтобы использовать ode15s или ode23t, вам нужны два указателя на функцию: одно представление большой матрицы системы ОДУ и другого представления вектора, содержащего правые стороны уравнений большой матрицы. Эти указатели на функцию являются эквивалентным представлением большой матрицы системы ОДУ где M (t, y (t)) y’ (t) = f (t, y (t)).

Найдите, что эти указатели на функцию при помощи massMatrixForm получают большую матрицу massM (M в уравнении) и правые стороны f.

[massM,f] = massMatrixForm(ODEs,vars)

massM =
[           -1,                     0,                     0,             0,             0]
[            0,                    -1,                     0,             0,             0]
[            0,                     0,                     0,             m,             0]
[            0,                     0,                     0,             0,             m]
[ -4*T(t)*x(t), 2*g*m*r - 4*T(t)*y(t), - 2*x(t)^2 - 2*y(t)^2, -4*m*r*Dxt(t), -4*m*r*Dyt(t)]
 
f =
               -Dxt(t)
               -Dyt(t)
         (T(t)*x(t))/r
 (T(t)*y(t) - g*m*r)/r
                     0

Уравнения в ODEs могут содержать символьные параметры, которые не заданы в векторе переменных vars. Найдите эти параметры при помощи setdiff на выводе symvar от ODEs и vars.

pODEs = symvar(ODEs);
pvars = symvar(vars);
extraParams = setdiff(pODEs, pvars)

extraParams =
[ g, m, r]

Дополнительными параметрами, которые необходимо задать, является массовый 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);

Шаг 3. Начальные условия для `ode15s` и `ode23t`

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

Во-первых, проверяйте переменные в vars.

vars

vars =
   x(t)
   y(t)
   T(t)
 Dxt(t)
 Dyt(t)

Здесь, 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)

Теперь создайте набор опции, который содержит большую матрицу M системы и векторный yp0 сопоставимых начальных значений для производных. Вы будете использовать этот набор опции при решении системы.

opt = odeset(opt, 'InitialSlope', yp0);

Шаг 4. Решите систему ОДУ с `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*cos(pi/6); -r*sin(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

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