Шаг 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. Уменьшение разностного порядка

Дифференциальный порядок системы дисковых полок является наивысшим дифференциальным порядком ее уравнений. Для решения дисковых полок с помощью 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. Уменьшить дифференциальный индекс с помощью reduceDAEToODE

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

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

ODEs = 
$$\left(\begin{array}{c}\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)\\ 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}\\ 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)-g\hspace{0.17em}m\hspace{0.17em}r}{r}\\ -\left(4\hspace{0.17em}T\left(t\right)\hspace{0.17em}y\left(t\right)-2\hspace{0.17em}g\hspace{0.17em}m\hspace{0.17em}r\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)-\left(2\hspace{0.17em}{x\left(t\right)}^2+2\hspace{0.17em}{y\left(t\right)}^2\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }T\left(t\right)-4\hspace{0.17em}T\left(t\right)\hspace{0.17em}x\left(t\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }x\left(t\right)-4\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxt}\left(t\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dxt}\left(t\right)-4\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dyt}\left(t\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)\end{array}\right)$$

constraints = 
$$\left(\begin{array}{c}-2\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}{\mathrm{Dxt}\left(t\right)}^2-2\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}{\mathrm{Dyt}\left(t\right)}^2-2\hspace{0.17em}T\left(t\right)\hspace{0.17em}{x\left(t\right)}^2-2\hspace{0.17em}T\left(t\right)\hspace{0.17em}{y\left(t\right)}^2+2\hspace{0.17em}g\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}y\left(t\right)\\ 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)\end{array}\right)$$

Шаг 4. ODE в дескрипторы функций для ode15 и 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 = 
$$\left(\begin{array}{ccccc}-1& 0& 0& 0& 0\\ 0& -1& 0& 0& 0\\ 0& 0& 0& m& 0\\ 0& 0& 0& 0& m\\ -4\hspace{0.17em}T\left(t\right)\hspace{0.17em}x\left(t\right)& 2\hspace{0.17em}g\hspace{0.17em}m\hspace{0.17em}r-4\hspace{0.17em}T\left(t\right)\hspace{0.17em}y\left(t\right)& -2\hspace{0.17em}{x\left(t\right)}^2-2\hspace{0.17em}{y\left(t\right)}^2& -4\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxt}\left(t\right)& -4\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dyt}\left(t\right)\end{array}\right)$$

f = 
$$\left(\begin{array}{c}-\mathrm{Dxt}\left(t\right)\\ -\mathrm{Dyt}\left(t\right)\\ \frac{T\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ \frac{T\left(t\right)\hspace{0.17em}y\left(t\right)-g\hspace{0.17em}m\hspace{0.17em}r}{r}\\ 0\end{array}\right)$$

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

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

extraParams = $$\left(\begin{array}{ccc}g& m& r\end{array}\right)$$

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

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

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

Сначала проверьте переменные в vars.

vars

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)$$

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

y0 = 5×1

    0.5000
   -0.8660
   -8.4957
         0
         0

yp0 = 5×1

         0
         0
         0
   -4.2479
   -2.4525

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

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

Шаг 6. Решение системы ОДУ с помощью 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