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

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

Создайте следующую систему двух дифференциальных алгебраических уравнений. Здесь символические функции 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