gusucode.com > symbolic工具箱matlab源码程序 > symbolic/@sym/findDecoupledBlocks.m

    function [eqsBlocks, varsBlocks] = findDecoupledBlocks(eqs, vars)
%FINDDECOUPLEDBLOCKS   Search for a block partitioning of 
%   equations to 'lower block triangular form'. For a square 
%   system of symbolic equations/expressions given in the 
%   vector eqs for the variables given in the vector vars,
%   the call
%
%   [eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs, vars)
%
%   serves for identifying subsets ('blocks') of equations 
%   that can be used to define subsets of the variables.
%   The cell arrays eqsBlocks and varsBlocks contain vectors
%   of indices, each of type double, that define the blocks. 
%   The i-th block of equations consists of the equations 
%   eqs(eqsBlocks{i}). It involves only the variables in 
%   vars(varsBlocks{1:i}). 
%
%   The i-th block is the set of equations determining the 
%   variables in vars(varsBlocks{i}). The parameters in 
%   vars(varsBlocks{i}),...,vars(varsBlocks{i-1}) were 
%   determined recursively by the previous blocks of equations.
%
%   In particular, when at least two blocks are found, the 
%   first block eqs(eqsBlocks{1}) defines a decoupled subset 
%   of equations containing only the subset of variables 
%   given by the first block of variables vars(varsBlock{1}).
%
%   After you solved 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 (which are known
%   at this time).
%
%   Thus, if a non-trivial block decomposition is possible,
%   you can split the solution process for a large system of 
%   equations involving a large number of variables into several
%   steps, where each step involves a smaller (sub-) system.
%
%   The number of blocks is given by length(eqsBlocks), which 
%   coincides with length(varsBlocks). If length(eqsBlocks) = 
%   length(varsBlocks) = 1, then no non-trivial block 
%   decomposition of the equations is possible.
%
%   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.
%
%   The implemented algorithm requires that for each variable
%   there is at least one 'matching' equation involving this 
%   variable that was not 'matched' to another variable. If 
%   the system does not satisfy this condition, then
%   FINDDECOUPLEDBLOCKS throws an error.
% 
%   In particular, length(eqs) must coincide with length(vars).
%  
%   Example:
%     Compute a lower triangular block decomposition of a 
%     symbolic system of differential algebraic equations
%     for the 'state variables' x1(t),x2(t),x3(t),x4(t). 
%     It contains constant symbolic parameters c1,c2,c3,c4 
%     and arbitrary symbolic functions f(t,x,y),g(t,x,y):
%
%     >> syms x1(t) x2(t) x3(t) x4(t) c1 c2 c3 c4 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)];
%        [eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs, vars);
%
%     The block decomposition returns a block of 2 equations 
%
%     >> eqs(eqsBlocks{1})
%
%     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))
%
%     for the 2 variables
%
%     >> vars(varsBlocks{1})
%
%     [x1(t),x3(t)]
%
%     After you solve this block for the variables x1(t),x3(t),
%     you can solve the next block of equations
%
%     >> eqs(eqsBlocks{2})
%
%     c2*diff(x1(t),t)+c1*diff(x3(t),t)==c4*g(t,x3(t),x4(t))
%
%     for the variable
%
%     >> vars(varsBlocks{2})
%
%     x4(t)
%
%     Once x4(t) is determined, the final block of equations
%
%     >> eqs(eqsBlocks{3})
%
%     x2(t)==f(t,x3(t),x4(t))
%
%     defines the variable
%
%     >> vars(varsBlocks{3})
%
%     x2(t)
%
%     Permuting the equations 
%     >> eqs = eqs([eqsBlocks{:}])
%
%     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))
%
%     and variables
%
%     >> vars = vars([varsBlocks{:}])
%
%     [x1(t),x3(t),x4(t),x2(t)]
%
%     produces a 'block lower triangular' system of equations
%     with three diagonal blocks of size 2x2,1x1, and 1x1:
%     >> incidenceMatrix(eqs, vars)
%
%     1     1     0     0
%     1     1     0     0
%     1     1     1     0
%     0     1     1     1
 
%   Copyright 2014 The MathWorks, Inc.

narginchk(2,2);
[eqs, vars] = checkDAEInput(eqs, vars);
blocks = feval(symengine, 'daetools::findDecoupledBlocks', eqs, vars);
blocks = num2cell(blocks);
eqsBlocks = cellfun(@double, num2cell(blocks{1}), 'Uniform', false);
varsBlocks = cellfun(@double, num2cell(blocks{2}), 'Uniform', false);
end