在matlab中求解矩阵DAE系统

Question

可以找到方程 here。如您所见，它是由 8 个标量方程组成的集合，接近 3 个矩阵方程。为了让 Matlab 知道方程是矩阵式的，我将可变时间相关向量函数声明为：

syms t p1(t) p2(t) p3(t)
p(t) = symfun([p1(t);p2(t);p3(t)], t);
p = formula(p(t));  % allows indexing for vector p
% same goes for w(t) and m(t)...

已知矩阵声明如下：

A = sym('A%d%d',[3 3]);
Fq = sym('Fq%d%d',[2 3]);
Im = diag(sym('Im%d%d',[1 3]));

系统现在可以根据guide:

建模

eqs = [diff(p) == A*w + Fq'*m,...
       diff(w) == -Im*p,...
       Fq*w == 0];
vars = [p; w; m];

此时，当我尝试减少索引（因为它等于 2）时，我收到以下错误：

[DAEs,DAEvars] = reduceDAEIndex(eqs,vars);

Error using sym/reduceDAEIndex (line 95)
Expecting as many equations as variables.

如果我们将所有变量都声明为标量，则不会出现错误：

syms A Im Fq real p(t) w(t) m(t)

引用 symfun documentation（提示部分）：

Symbolic functions are always scalars, therefore, you cannot index into a function.

然而，我很难相信不可能以矩阵方式求解这些方程。显然，可以将其扩展为 8 个标量方程，但这里涉及的多体系统非常简单，目的是能够求解复杂的 - 因此问题：是否可以在 Matlab 中求解矩阵 DAE，如果可以- 必须解决什么问题才能使其正常工作？

Ps。 Matlab DAE 求解器还有另一个问题：我的模型的输入变量（已知系数函数）是时变的。就 example is concerned, they are constant in all domain, however for my problem they change in time. This problem has been brought out here 而言。如有解决办法，请指教。

Answer 1

最后，我设法找到了这个问题的正确语法。我错误地将矩阵变量（例如 A、Fq）视为单个实体。下面我展示了利用矩阵方法解决这个特定 DAE 的代码：

% Define symbolic variables.
q = sym('q%d',[3 1]);        % state variables
a = sym('a'); k = sym('k');  % constant parameters
t = sym('t','real');         % independent variable

% Define system variables and group them in vectors:
p1(t) = sym('p1(t)'); p2(t) = sym('p2(t)'); p3(t) = sym('p3(t)');
w1(t) = sym('w1(t)'); w2(t) = sym('w2(t)'); w3(t) = sym('w3(t)');
m1(t) = sym('m1(t)'); m2(t) = sym('m2(t)');
pvect = [p1(t); p2(t); p3(t)];
wvect = [w1(t); w2(t); w3(t)];
mvect = [m1(t); m2(t)];

% Define matrices:
mass = diag(sym('ms%d',[1 3]));
Fq = [0 -1 a;
      0  0 1];
A = [1  0  0;
     0  1  a;
     0  a -q(1)*a] * k;

% Define sets of equations and aggregate them into one set: 
set1 = diff(pvect,t) == A*wvect + Fq'*mvect;
set2 = mass*diff(wvect,t) == -pvect;
set3 = Fq*wvect == 0;
eqs = [set1; set2; set3];
% Close all system variables in one vector:
vars = [pvect; wvect; mvect];

% Reduce index of the system and remove redundnat equations:
[DAEs,DAEvars] = reduceDAEIndex(eqs,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
[M,F] = massMatrixForm(DAEs,DAEvars);

我们收到两个变量 p1(t) 和 w1(t) 的非常简单的 2x2 ODE。请记住，在减少冗余后，我们摆脱了状态向量 q 中的所有元素。这意味着所有左变量（k 和 mass(1,1)）都与时间无关。如果系统中的某些变量存在时间依赖性，那么这个案子就更难解决了。

% Replace symbolic variables with numeric ones:
M = odeFunction(M, DAEvars,mass(1,1));
F = odeFunction(F, DAEvars, k);
k = 2000; numericMass = 4;
F = @(t, Y) F(t, Y, k);
M = @(t, Y) M(t, Y, numericMass);

% set the solver:
opt = odeset('Mass', M); % Mass matrix of the system
TIME = [1; 0];           % Time boundaries of the simulation (backwards in time)
y0 = [1 0]';             % Initial conditions for left variables p1(t) and w1(t)

% Call the solver
[T, solution] = ode15s(F, TIME, y0, opt);
% Plot results
plot(T,solution(:,1),T,solution(:,2))

在matlab中求解矩阵DAE系统

Solve matrix DAE system in matlab

matlab

matrix

symbolic-math

ode

differential-equations