Main Content

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`