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

Этот рабочий процесс является альтернативным рабочим процессом к решению полулинейных дифференциальных алгебраических уравнений (ДАУ), используемые только если reduceDAEIndex отказавший в стандартном рабочем процессе с предупреждающим сообщением: The index of the reduced DAEs is larger than 1. [daetools::reduceDAEIndex]. Для стандартного рабочего процесса смотрите, Решают Дифференциальные Алгебраические уравнения (ДАУ).

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

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

Шаг 1. Задайте уравнения и переменные

Система уравнений ДАУ:

md2xdt2=T(t)x(t)rmd2ydt2=T(t)y(t)r-mgx(t)2+y(t)2=r2

Задайте независимые переменные и переменные состояния при помощи 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 = 

(mt Dxt(t)-T(t)x(t)rgm+mt Dyt(t)-T(t)y(t)r-r2+x(t)2+y(t)2Dxt(t)-t x(t)Dyt(t)-t y(t))[m*diff (Dxt (t), t) - (T (t) *x (t))/r; g*m + m*diff (Dyt (t), t) - (T (t) *y (t))/r; - r^2 + x (t) ^2 + y (t) ^2; Dxt (t) - diff (x (t), t); Dyt (t) - diff (y (t), t)]

vars = 

(x(t)y(t)T(t)Dxt(t)Dyt(t))xt ; yt ; T (t); Dxt (t); Dyt (t)]

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

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

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

(Dxt(t)-t x(t)Dyt(t)-t y(t)mt Dxt(t)-T(t)x(t)rmt Dyt(t)-T(t)y(t)-gmrr-4T(t)y(t)-2gmrt y(t)-2x(t)2+2y(t)2t T(t)-4T(t)x(t)t x(t)-4mrDxt(t)t Dxt(t)-4mrDyt(t)t Dyt(t))[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 (т) *y (t) - 2*g*m*r) *diff (y (t), t) - (2*x (т) ^2 + 2*y (т) ^2) *diff (T (t), t) - 4*T (т) *x (t) *diff (x (t), t) - 4*m*r*Dxt (т) *diff (Dxt (t), t) - 4*m*r*Dyt (т) *diff (Dyt (t), t)]

constraints = 

(-2mrDxt(t)2-2mrDyt(t)2-2T(t)x(t)2-2T(t)y(t)2+2gmry(t)r2-x(t)2-y(t)22Dxt(t)x(t)+2Dyt(t)y(t))[-2*m*r*Dxt (т) ^2 - 2*m*r*Dyt (т) ^2 - 2*T (т) *x (t) ^2 - 2*T (т) *y (t) ^2 + 2*g*m*r*y (т); r^2 - x (t) ^2 - y (t) ^2; 2*Dxt (т) *x (t) + 2*Dyt (т) *y (t)]

Шаг 4. ОДУ к Указателям на функции для 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 = 

(-100000-1000000m00000m-4T(t)x(t)2gmr-4T(t)y(t)-2x(t)2-2y(t)2-4mrDxt(t)-4mrDyt(t))[-sym (1), sym (0), sym (0), sym (0), sym (0); sym (0),-sym (1), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), m, sym (0); sym (0), sym (0), sym (0), sym (0), m;-4*T (t) *x (t), 2*g*m*r - 4*T (т) *y (t), - 2*x (т) ^2 - 2*y (т) ^2,-4*m*r*Dxt (t),-4*m*r*Dyt (t)]

f = 

(-Dxt(t)-Dyt(t)T(t)x(t)rT(t)y(t)-gmrr0)[-Dxt (t);-Dyt (t); (T (t) *x (t))/r; (T (t) *y (t) - g*m*r)/r; sym (0)]

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

pODEs = symvar(ODEs);
pvars = symvar(vars);
extraParams = setdiff(pODEs, pvars)
extraParams = (gmr)[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);

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

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

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

vars
vars = 

(x(t)y(t)T(t)Dxt(t)Dyt(t))xt ; yt ; 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)
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

Решите систему, объединяющуюся по отрезку времени 0t≤ 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

Figure contains an axes. The axes contains 3 objects of type line. These objects represent x(t), y(t), T(t).

Решите систему для различных значений параметров путем устанавливания нового значения и регенерации указателя на функцию и начальных условий.

Установите 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

Figure contains an axes. The axes contains 3 objects of type line. These objects represent x(t), y(t), T(t).

Смотрите также

| | | | | | | | | |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте