Блок нижнего треугольного разложения системы ДАУ

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

Создайте следующую систему четырех дифференциальных алгебраических уравнений. Здесь символическая функция вызывает x1(t), x2(t), x3(t), и x4(t) представляет переменные состояния системы. Система также содержит символические параметры c1, c2, c3, c4, и функции f(t,x,y) и g(t,x,y).

syms x1(t) x2(t) x3(t) x4(t)
syms c1 c2 c3 c4
syms f(t,x,y) g(t,x,y)

eqs = [c1*diff(x1(t),t)+c2*diff(x3(t),t)==c3*f(t,x1(t),x3(t));...
       c2*diff(x1(t),t)+c1*diff(x3(t),t)==c4*g(t,x3(t),x4(t));...
       x1(t)==g(t,x1(t),x3(t));...
       x2(t)==f(t,x3(t),x4(t))];

vars = [x1(t), x2(t), x3(t), x4(t)];

Использовать findDecoupledBlocks для поиска блока структуры системы.

[eqsBlocks, varsBlocks] = findDecoupledBlocks(eqs, vars)

eqsBlocks =
  1×3 cell array
    {1×2 double}    {[2]}    {[4]}
varsBlocks =
  1×3 cell array
    {1×2 double}    {[4]}    {[2]}

Первый блок содержит два уравнения в двух переменных.

eqs(eqsBlocks{1})

ans =
 c1*diff(x1(t), t) + c2*diff(x3(t), t) == c3*f(t, x1(t), x3(t))
                                    x1(t) == g(t, x1(t), x3(t))

vars(varsBlocks{1})

ans =
[ x1(t), x3(t)]

После того, как вы решите этот блок для переменных x1(t), x3(t)можно решить следующий блок уравнений. Этот блок состоит из одного уравнения.

eqs(eqsBlocks{2})

ans =
c2*diff(x1(t), t) + c1*diff(x3(t), t) == c4*g(t, x3(t), x4(t))

Блок включает одну переменную.

vars(varsBlocks{2})

ans =
x4(t)

После того, как вы решите уравнение из блока 2 для переменной x4(t), оставшийся блок уравнений, eqs(eqsBlocks{3}), задает оставшуюся переменную, vars(varsBlocks{3}).

eqs(eqsBlocks{3})
vars(varsBlocks{3})

ans =
x2(t) == f(t, x3(t), x4(t))
 
ans =
x2(t)

Найдите сочетания, которые преобразуют систему в блок нижнюю треугольную форму.

eqsPerm = [eqsBlocks{:}]
varsPerm = [varsBlocks{:}]

eqsPerm =
     1     3     2     4

varsPerm =
     1     3     4     2

Преобразуйте систему в блочную нижнюю треугольную систему уравнений.

eqs = eqs(eqsPerm)
vars = vars(varsPerm)

eqs =
 c1*diff(x1(t), t) + c2*diff(x3(t), t) == c3*f(t, x1(t), x3(t))
                                    x1(t) == g(t, x1(t), x3(t))
 c2*diff(x1(t), t) + c1*diff(x3(t), t) == c4*g(t, x3(t), x4(t))
                                    x2(t) == f(t, x3(t), x4(t))

vars =
[ x1(t), x3(t), x4(t), x2(t)]

Найдите матрицу инцидентности получившейся системы. Матрица инцидентности показывает, что система перестановочных уравнений имеет три диагональных блока размера 2-by- 2, 1-by- 1, и 1-by- 1.

incidenceMatrix(eqs, vars)

ans =
     1     1     0     0
     1     1     0     0
     1     1     1     0
     0     1     1     1

Разложение BLT и решение линейной системы

Найдите блоки уравнений в линейной алгебраической системе, а затем решите систему, последовательно решая каждый блок уравнений, начиная с первого.

Создайте следующую систему линейных алгебраических уравнений.

syms x1 x2 x3 x4 x5 x6 c1 c2 c3

eqs = [c1*x1 + x3 + x5 == c1 + c2 + 1;...
       x1 + x3 + x4 + 2*x6 == 4 + c2;...
       x1 + 2*x3 + c3*x5 == 1 + 2*c2 + c3;...
       x2 + x3 + x4 + x5 == 2 + c2;...
       x1 - c2*x3 + x5 == 2 - c2^2;...
       x1 - x3 + x4 - x6 == 1 - c2];

vars = [x1, x2, x3, x4, x5, x6];

Использовать findDecoupledBlocks для преобразования системы в более низкую треугольную форму. Для этой системы, findDecoupledBlocks определяет три блока уравнений и соответствующих переменных.

[eqsBlocks, varsBlocks] = findDecoupledBlocks(eqs, vars)

eqsBlocks =
  1×3 cell array
    {1×3 double}    {1×2 double}    {[4]}
varsBlocks =
  1×3 cell array
    {1×3 double}    {1×2 double}    {[2]}

Идентифицируйте переменные в первом блоке. Этот блок состоит из трех уравнений в трех переменных.

vars(varsBlocks{1})

ans =
[ x1, x3, x5]

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

[x1,x3,x5] = solve(eqs(eqsBlocks{1}), vars(varsBlocks{1}))

x1 =
1
 
x3 =
c2
 
x5 =
1

Идентифицируйте переменные во втором блоке. Этот блок состоит из двух уравнений в двух переменных.

vars(varsBlocks{2})

ans =
[ x4, x6]

Решите этот блок уравнений, присваивая решения соответствующим переменным.

[x4, x6] = solve(eqs(eqsBlocks{2}), vars(varsBlocks{2}))

x4 =
x3/3 - x1 - c2/3 + 2
 
x6 =
(2*c2)/3 - (2*x3)/3 + 1

Использовать subs чтобы вычислить результат для уже известных значений переменных x1, x3, и x5.

x4 = subs(x4)
x6 = subs(x6)

x4 =
1
 
x6 =
1

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

vars(varsBlocks{3})

ans =
x2

Решите это уравнение, присвоив решение x2.

x2 = solve(eqs(eqsBlocks{3}), vars(varsBlocks{3}))

x2 =
c2 - x3 - x4 - x5 + 2

Использовать subs чтобы вычислить результат для уже известных значений всех других переменных этой системы.

x2 = subs(x2)

x2 =
0

Также можно переписать этот пример с помощью for-цикл. Этот подход позволяет вам использовать пример для больших систем уравнений.

syms x1 x2 x3 x4 x5 x6 c1 c2 c3

eqs = [c1*x1 + x3 + x5 == c1 + c2 + 1;...
       x1 + x3 + x4 + 2*x6 == 4 + c2;...
       x1 + 2*x3 + c3*x5 == 1 + 2*c2 + c3;...
       x2 + x3 + x4 + x5 == 2 + c2;...
       x1 - c2*x3 + x5 == 2 - c2^2
       x1 - x3 + x4 - x6 == 1 - c2];

vars = [x1, x2, x3, x4, x5, x6];

[eqsBlocks, varsBlocks] = findDecoupledBlocks(eqs, vars);

vars_sol = vars;

for i = 1:numel(eqsBlocks)
  sol = solve(eqs(eqsBlocks{i}), vars(varsBlocks{i}));
  vars_sol_per_block = subs(vars(varsBlocks{i}), sol);
  for k=1:i-1
    vars_sol_per_block = subs(vars_sol_per_block, vars(varsBlocks{k}),...
                         vars_sol(varsBlocks{k}));
  end
  vars_sol(varsBlocks{i}) = vars_sol_per_block
end

vars_sol =
[ 1, x2, c2, x4, 1, x6]
 
vars_sol =
[ 1, x2, c2, 1, 1, 1]
 
vars_sol =
[ 1, 0, c2, 1, 1, 1]