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

Этот рабочий процесс является альтернативным рабочим процессом к решению полулинейных ДАУ, используемых только если 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 переменных- 1 вектор. Поэтому предположениями для начальных значений переменных и их производных должен также быть 5- 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*sin(pi/6); -r*cos(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

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