Блок нижней треугольной декомпозиции системы дисковых полок

Вычислите блок нижней треугольной декомпозиции (БЛТ декомпозиции) символической системы дифференциальных алгебраических уравнений (DAE).

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

incidenceMatrix(eqs, vars)

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

Декомпозиция БЛТ и решение линейной системы

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

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

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]