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

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

Создайте следующую систему четырех дифференциальных алгебраических уравнений. Здесь, символьные вызовы функции x1(t), x2(t), x3(t), и x4(t) представляйте переменные состояния системы. Система также содержит символьные параметры c1C2 , 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

Разложение 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]