Search for decoupled blocks in systems of equations

`[`

identifies
subsets (blocks) of equations that can be used to define subsets of
variables. The number of variables `eqsBlocks`

,`varsBlocks`

]
= findDecoupledBlocks(`eqs`

,`vars`

)`vars`

must
coincide with the number of equations `eqs`

.

The *i*th 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.

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

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]

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.

`daeFunction`

| `decic`

| `diag`

| `incidenceMatrix`

| `isLowIndexDAE`

| `massMatrixForm`

| `odeFunction`

| `reduceDAEIndex`

| `reduceDAEToODE`

| `reduceDifferentialOrder`

| `reduceRedundancies`

| `tril`

| `triu`