Преобразуйте систему ДАУ в неявную систему ОДУ

Преобразуйте систему дифференциальных алгебраических уравнений (ДАУ) к системе неявных обыкновенных дифференциальных уравнений (ОДУ).

Создайте следующую систему двух дифференциальных алгебраических уравнений. Здесь, символьные функции x(t), y(t) и z(t) представляют переменные состояния системы. Задайте уравнения и переменные как два символьных вектора: уравнения как вектор символьных уравнений и переменные как вектор символьных вызовов функции.

syms x(t) y(t) z(t)
eqs = [diff(x,t)+x*diff(y,t) == y,...
       x*diff(x, t)+x^2*diff(y) == sin(x),...
       x^2 + y^2 == t*z];
vars = [x(t), y(t), z(t)];

Используйте reduceDAEToODE, чтобы переписать систему так, чтобы дифференциальным индексом был 0.

newEqs = reduceDAEToODE(eqs, vars)

newEqs =
                            x(t)*diff(y(t), t) - y(t) + diff(x(t), t)
                diff(x(t), t)*(cos(x(t)) - y(t)) - x(t)*diff(y(t), t)
 z(t) - 2*x(t)*diff(x(t), t) - 2*y(t)*diff(y(t), t) + t*diff(z(t), t)

Уменьшайте систему и возвратите больше деталей

Проверяйте, имеет ли следующая система ДАУ минимум (0 или 1) или высокий (>1) дифференциальный индекс. Если индекс выше, чем 1, сначала попытайтесь уменьшать индекс при помощи reduceDAEIndex и затем при помощи reduceDAEToODE.

Создайте систему дифференциальных алгебраических уравнений. Здесь, функции x1(t), x2(t) и x3(t) представляют переменные состояния системы. Система также содержит функции q1(t), q2(t) и q3(t). Эти функции не представляют переменные состояния. Задайте уравнения и переменные как два символьных вектора: уравнения как вектор символьных уравнений и переменные как вектор символьных вызовов функции.

syms x1(t) x2(t) x3(t) q1(t) q2(t) q3(t)
eqs = [diff(x2) == q1 - x1,
       diff(x3) == q2 - 2*x2 - t*(q1-x1),
       q3 - t*x2 - x3];
vars = [x1(t), x2(t), x3(t)];

Используйте isLowIndexDAE, чтобы проверять дифференциальный индекс системы. Для этой системы isLowIndexDAE возвращает 0 (false). Это означает, что дифференциальным индексом системы является 2 или выше.

isLowIndexDAE(eqs, vars)

ans =
  logical
     0

Используйте reduceDAEIndex в качестве своей первой попытки переписать систему так, чтобы дифференциальным индексом был 1. Для этой системы reduceDAEIndex выдает предупреждение, потому что это не может уменьшать дифференциальный индекс системы к 0 или 1.

[newEqs, newVars] = reduceDAEIndex(eqs, vars)

Warning: Index of reduced DAEs is larger than 1.
 
newEqs =
                      x1(t) - q1(t) + diff(x2(t), t)
       Dx3t(t) - q2(t) + 2*x2(t) + t*(q1(t) - x1(t))
                             q3(t) - x3(t) - t*x2(t)
 diff(q3(t), t) - x2(t) - t*diff(x2(t), t) - Dx3t(t)
 
newVars =
   x1(t)
   x2(t)
   x3(t)
 Dx3t(t)

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

[newEqs, constraintEqs] = reduceDAEToODE(eqs, vars)

newEqs =
                                             x1(t) - q1(t) + diff(x2(t), t)
                       2*x2(t) - q2(t) + t*q1(t) - t*x1(t) + diff(x3(t), t)
 diff(x1(t), t) - diff(q1(t), t) + diff(q2(t), t, t) - diff(q3(t), t, t, t)

constraintEqs =
 x1(t) - q1(t) + diff(q2(t), t) - diff(q3(t), t, t)
                            x3(t) - q3(t) + t*x2(t)
                     x2(t) - q2(t) + diff(q3(t), t)

Используйте синтаксис с тремя выходными аргументами, чтобы возвратить новые уравнения, уравнения ограничений и дифференциальный индекс исходной системы, eqs.

[newEqs, constraintEqs, oldIndex] = reduceDAEToODE(eqs, vars)

newEqs =
                                             x1(t) - q1(t) + diff(x2(t), t)
                       2*x2(t) - q2(t) + t*q1(t) - t*x1(t) + diff(x3(t), t)
 diff(x1(t), t) - diff(q1(t), t) + diff(q2(t), t, t) - diff(q3(t), t, t, t)

constraintEqs =
 x1(t) - q1(t) + diff(q2(t), t) - diff(q3(t), t, t)
                            x3(t) - q3(t) + t*x2(t)
                     x2(t) - q2(t) + diff(q3(t), t)

oldIndex =
     3