Résoudre le système DAE matriciel dans matlab

Question

Les équations peuvent être trouvées here. Comme vous pouvez le voir, il s'agit d'un ensemble de 8 équations scalaires fermées à 3 matrices. Afin de laisser Matlab savoir que les équations sont matrice - sage, je déclare les fonctions de temps variable de vecteur dépendant comme:

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)... 

matrices connues sont déclarées comme suit:

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

Le système est maintenant prêt à être modélisé selon guide:

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

A ce moment, lorsque je tente de réduire l'indice (puisqu'il est égal à 2), je reçois l'erreur suivante:

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

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

---

L'erreur ne se poserait pas si nous avions déclaré toutes les variables comme scalaires:

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

Citant symfundocumentation (section Conseils):

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

Cependant, il est difficile pour moi de croire qu'il n'est pas possible de résoudre ces équations matricielles. Évidemment, on peut l'étendre à 8 équations scalaires, mais le système multi-corps concerné ici est très simple et le but est de pouvoir résoudre des équations complexes - d'où la question: est-il possible de résoudre la matrice DAE dans Matlab, et si oui - Qu'est-ce qui doit être réparé pour que cela fonctionne?

---

Ps. J'ai un autre problème avec le solveur Matlab DAE: les variables d'entrée (fonctions de coefficients connues) pour mon modèle sont des variantes temporelles. En ce qui concerne example, ils sont constants dans tous les domaines, mais pour mon problème, ils changent dans le temps. Ce problème a été mis en évidence here. Je vous serais reconnaissant si vous en parliez, si vous aviez une solution.

Answer 1

Enfin, j'ai réussi à trouver la syntaxe correcte pour ce problème. J'ai fait une erreur de traitement des variables matricielles (telles que A, Fq) en tant qu'entité unique. Ci-dessous je vous présente le code qui utilise l'approche matricielle et permet de résoudre ce DAE particulier:

% 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); 

Nous recevons très simple 2x2 ODE pour deux variables p1(t) et w1(t). Gardez à l'esprit qu'après avoir réduit les redondances, nous nous sommes débarrassés de tous les éléments du vecteur d'état q. Cela signifie que toutes les variables de gauche (k et mass(1,1)) ne dépendent pas du temps. S'il y avait eu une dépendance temporelle de certaines variables dans le système, le cas aurait été beaucoup plus difficile à résoudre.

% 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))