findDecoupledBlocks

Search for decoupled blocks in systems of equations

Syntax

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

Description

[eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs,vars) identifies subsets (blocks) of equations that can be used to define subsets of variables. The number of variables vars must coincide with the number of equations eqs.

The ith block is the set of equations determining the variables in vars(varsBlocks{i}). The variables in vars([varsBlocks{1},…,varsBlocks{i-1}]) are determined recursively by the previous blocks of equations. After you solve the first block of equations for the first block of variables, the second block of equations, given by eqs(eqsBlocks{2}), defines a decoupled subset of equations containing only the subset of variables given by the second block of variables, vars(varsBlock{2}), plus the variables from the first block (these variables are known at this time). Thus, if a nontrivial block decomposition is possible, you can split the solution process for a large system of equations involving many variables into several steps, where each step involves a smaller subsystem.

The number of blocks length(eqsBlocks) coincides with length(varsBlocks). If length(eqsBlocks) = length(varsBlocks) = 1, then a nontrivial block decomposition of the equations is not possible.

Examples

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]

Input Arguments

eqs — System of equations

vector of symbolic equations | vector of symbolic expressions

System of equations, specified as a vector of symbolic equations or expressions.

vars — Variables

vector of symbolic variables | vector of symbolic functions | vector of symbolic function calls

Variables, specified as a vector of symbolic variables, functions, or function calls, such as x(t).

Example: [x(t),y(t)] or [x(t);y(t)]

Output Arguments

eqsBlocks — Indices defining blocks of equations

cell array

Indices defining blocks of equations, returned as a cell array. Each block of indices is a row vector of double-precision integer numbers. The ith block of equations consists of the equations eqs(eqsBlocks{i}) and involves only the variables in vars(varsBlocks{1:i}).

varsBlocks — Indices defining blocks of variables

cell array

Indices defining blocks of variables, returned as a cell array. Each block of indices is a row vector of double-precision integer numbers. The ith block of equations consists of the equations eqs(eqsBlocks{i}) and involves only the variables in vars(varsBlocks{1:i}).

Tips

• The implemented algorithm requires that for each variable in vars there must be at least one matching equation in eqs involving this variable. The same equation cannot also be matched to another variable. If the system does not satisfy this condition, then findDecoupledBlocks throws an error. In particular, findDecoupledBlocks requires that length(eqs) = length(vars).

• Applying the permutations e = [eqsBlocks{:}] to the vector eqs and v = [varsBlocks{:}] to the vector vars produces an incidence matrix incidenceMatrix(eqs(e), vars(v)) that has a block lower triangular sparsity pattern.