Решение дифференциальных алгебраических уравнений (дисковых полок)

Открыть сценарий в реальном времени

В этом примере показано, как решать дифференциальные алгебраические уравнения (дисковые полки) с помощью Toolbox™ MATLAB ® и Symbolic Math.

Дифференциальные алгебраические уравнения, включающие функции или переменные состояния, $x (t) =_{} [x1 (t), ._{.} .,$ xn (t)] имеют вид

$F (t, x (t_{}^{)}, x˙ ($ t)) = 0

где $t$ - независимая переменная. Число уравнений $F =_{[} F1, . . ._{,}$ Fn] должно соответствовать числу $переменных$ состояния x (t) = [x1 (t),..., xn (t)].

Поскольку большинство систем дисковых полок не подходят для прямого ввода в решатели MATLAB ®, такие какode15iсначала преобразуйте их в подходящую форму с помощью функции символьных математических Toolbox™. Эта функциональность уменьшает дифференциальный индекс (количество дифференциаций, необходимых для уменьшения системы до ОДУ) дисковых полок до 1 или 0, а затем преобразует систему дисковых полок в числовые функциональные дескрипторы, подходящие для решателей MATLAB ®. Затем используйте решатели MATLAB ®, такие как ode15i, ode15s, или ode23t, для решения дисковых полок.

Выполните следующие действия для решения проблемы системы дисковых полок.

Шаг 1: Определение уравнений и переменных

На следующем рисунке показан рабочий процесс дисковой полки путем решения дисковых полок для маятника.

Переменные состояния:

Горизонтальное положение маятника $x (t$ )
Вертикальное положение маятника $y (t$ )
Сила, предотвращающая отлет маятника $от T (t$ )

Переменные:

Маятниковая масса $м$
Длина маятника $r$
Гравитационная постоянная $g$

Система уравнений дисковых полок:

$\begin{array}{l} \frac{^{}}{^{md2xdt2}} = T \frac{(t)}{x} \\ \frac{(^{t})}{^{}} \frac{rmd2ydt2 =}{T} (t \\ ) y^{(} t)^{} {r-mgx}^{} \end{array}$ (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: Уменьшение дифференциального порядка

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

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ T (t) \\ Dxt (t) \\ Dyt (t) \end{array})$

Шаг 3: Проверка и уменьшение дифференциального индекса

3.1 Проверить дифференциальный индекс системы

Проверьте дифференциальный индекс системы дисковой полки с помощью isLowIndexDAE. Если индекс равен 0 или 1, то isLowIndexDAE возвращает логический 1 (true) и можно пропустить шаг 3.2 и перейти к шагу 4. Преобразование систем дисковых полок в дескрипторы функций MATLAB. Здесь, isLowIndexDAE возвращает логический 0 (false), что означает, что дифференциальный индекс больше, чем 1 и должны быть сокращены.

isLowIndexDAE(eqns,vars)

ans = logical
   0

3.2 Уменьшить дифференциальный индекс с помощью reduceDAEIndex

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

Уменьшите дифференциальный индекс дисковых полок, описанный в eqns и vars.

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)

DAEs = 
 $(\begin{array}{c} m Dxtt (t) - \frac{T (t) x (t)}{r} \\ g m + m Dytt (t) - \frac{T (t) y (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ Dxt (t) - {Dxt}_{1} (t) \\ Dyt (t) - {Dyt}_{1} (t) \\ 2 {Dxt}_{1} (t) x (t) + 2 {Dyt}_{1} (t) y (t) \\ 2 y (t) \frac{\partial}{\partial t} {Dyt}_{1} (t) + 2 {{Dxt}_{1} (t)}^{2} + 2 {{Dyt}_{1} (t)}^{2} + 2 Dxt1t (t) x (t) \\ Dxtt (t) - Dxt1t (t) \\ Dytt (t) - \frac{\partial}{\partial t} {Dyt}_{1} (t) \\ {Dyt}_{1} (t) - \frac{\partial}{\partial t} y (t) \end{array})$

DAEvars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ T (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \\ {Dxt}_{1} (t) \\ {Dyt}_{1} (t) \\ Dxt1t (t) \end{array})$

Если reduceDAEIndex выдает ошибку или предупреждение, используйте альтернативный рабочий процесс, описанный в разделе Решение системы полулинейных дисковых полок.

Часто, reduceDAEIndex вводит избыточные уравнения и переменные, которые можно исключить. Устранение избыточных уравнений и переменных с помощью reduceRedundancies.

[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)

DAEs = 
 $(\begin{array}{c} - \frac{T (t) x (t) - m r Dxtt (t)}{r} \\ \frac{g m r - T (t) y (t) + m r Dytt (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ 2 Dxt (t) x (t) + 2 Dyt (t) y (t) \\ 2 {Dxt (t)}^{2} + 2 {Dyt (t)}^{2} + 2 Dxtt (t) x (t) + 2 Dytt (t) y (t) \\ Dytt (t) - \frac{\partial}{\partial t} Dyt (t) \\ Dyt (t) - \frac{\partial}{\partial t} y (t) \end{array})$

DAEvars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ T (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \end{array})$

Проверьте дифференциальный индекс новой системы. Сейчас, isLowIndexDAE возвращает логический 1 (true), что означает, что дифференциальный индекс системы равен 0 или 1.

isLowIndexDAE(DAEs,DAEvars)

ans = logical
   1

Шаг 4: Преобразование систем дисковых полок в дескрипторы функций MATLAB

На этом шаге создаются дескрипторы функций для решателя ODE MATLAB ®ode15i, который является решателем общего назначения. Использование специализированных решателей массовых матриц, таких как ode15s и ode23t, см. разделы Решение дисковых полок с помощью решателей массовых матриц и Выбор решателя ОДУ.

reduceDAEIndex выводит вектор уравнений в DAEs и вектор переменных в DAEvars. Использовать ode15i, необходим функциональный дескриптор, описывающий систему дисковых полок.

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

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars)

extraParams =  $(\begin{array}{c} g & m & r \end{array})$

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

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

DAEvars

DAEvars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ T (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \end{array})$

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

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

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;
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

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

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

Решение дифференциальных алгебраических уравнений (дисковых полок)

Шаг 1: Определение уравнений и переменных

Шаг 2: Уменьшение дифференциального порядка

Шаг 3: Проверка и уменьшение дифференциального индекса

Шаг 4: Преобразование систем дисковых полок в дескрипторы функций MATLAB

Шаг 5: Поиск начальных условий для решателей

Шаг 6: Решение дисковых полок с помощью `ode15i`

Связанные темы

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

Поддержка

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

Решение дифференциальных алгебраических уравнений (дисковых полок)

Шаг 1: Определение уравнений и переменных

Шаг 2: Уменьшение дифференциального порядка

Шаг 3: Проверка и уменьшение дифференциального индекса

Шаг 4: Преобразование систем дисковых полок в дескрипторы функций MATLAB

Шаг 5: Поиск начальных условий для решателей

Шаг 6: Решение дисковых полок с помощью ode15i

Связанные темы

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

Поддержка

Шаг 6: Решение дисковых полок с помощью `ode15i`