这些方程式可以在这里找到。如您所见，它由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]));

现在可以根据指南对系统进行建模：

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文档(提示部分)：

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

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

附言Matlab DAE求解器还有另一个问题：模型的输入变量(已知系数函数)是时变的。就示例而言，它们在所有域中都是恒定的，但是对于我的问题，它们随时间变化。这个问题已经在这里提出了。如果您有任何解决方案，请参考。

最后，我设法找到了解决此问题的正确语法。 我将矩阵变量(例如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))