Block Lower Triangular Decomposition of DAE System

Compute a block lower triangular decomposition (BLT decomposition) of a symbolic system of differential algebraic equations (DAEs).

Create the following system of four differential algebraic equations. Here, the symbolic function calls x1(t), x2(t), x3(t), and x4(t) represent the state variables of the system. The system also contains symbolic parameters c1, c2, c3, c4, and functions f(t,x,y) and 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)];

Use findDecoupledBlocks to find the block structure of the system.

[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]}

The first block contains two equations in two variables.

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)]

After you solve this block for the variables x1(t), x3(t), you can solve the next block of equations. This block consists of one equation.

eqs(eqsBlocks{2})

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

The block involves one variable.

vars(varsBlocks{2})

ans =
x4(t)

After you solve the equation from block 2 for the variable x4(t), the remaining block of equations, eqs(eqsBlocks{3}), defines the remaining variable, vars(varsBlocks{3}).

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

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

Find the permutations that convert the system to a block lower triangular form.

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

eqsPerm =
     1     3     2     4

varsPerm =
     1     3     4     2

Convert the system to a block lower triangular system of equations.

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)]

Find the incidence matrix of the resulting system. The incidence matrix shows that the system of permuted equations has three diagonal blocks of size 2-by-2, 1-by-1, and 1-by-1.

incidenceMatrix(eqs, vars)

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

BLT Decomposition and Solution of Linear System

Find blocks of equations in a linear algebraic system, and then solve the system by sequentially solving each block of equations starting from the first one.

Create the following system of linear algebraic equations.

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];

Use findDecoupledBlocks to convert the system to a lower triangular form. For this system, findDecoupledBlocks identifies three blocks of equations and corresponding variables.

[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]}

Identify the variables in the first block. This block consists of three equations in three variables.

vars(varsBlocks{1})

ans =
[ x1, x3, x5]

Solve the first block of equations for the first block of variables assigning the solutions to the corresponding variables.

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

x1 =
1
 
x3 =
c2
 
x5 =
1

Identify the variables in the second block. This block consists of two equations in two variables.

vars(varsBlocks{2})

ans =
[ x4, x6]

Solve this block of equations assigning the solutions to the corresponding variables.

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

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

Use subs to evaluate the result for the already known values of variables x1, x3, and x5.

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

x4 =
1
 
x6 =
1

Identify the variables in the third block. This block consists of one equation in one variable.

vars(varsBlocks{3})

ans =
x2

Solve this equation assigning the solution to x2.

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

x2 =
c2 - x3 - x4 - x5 + 2

Use subs to evaluate the result for the already known values of all other variables of this system.

x2 = subs(x2)

x2 =
0

Alternatively, you can rewrite this example using the for-loop. This approach lets you use the example for larger systems of equations.

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]