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

Вычислите блок нижнее треугольное разложение (разложение 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]