DAE Error using mupadengin​e/feval_in​ternal Expecting as many equations as variables.

I have a set of 15 equations:

eqns = [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10, eqn11, eqn12, eqn13, eqn14, eqn15];

I also have 15 variables, which are used in the equations:

p_vars = [p1(t); p2(t); p3(t); p4(t); p5(t); p6(t)];
q_vars = [q1(t); q2(t); q3(t); q4(t); q5(t); q6(t); q7(t); q8(t); q9(t)];
vars = [q_vars; p_vars]

Using isLowIndexDAE returned 0 so I decided to run reduceDAEIndex:

[newEqns, newVars, ~, oldIndex] = reduceDAEIndex(h_eqns, vars)

However it gives me an error, saying "Expecting as many equations as variables". This confuses me because it seems that there obviously are as many equations as variables.

Does anyone know why this is happening?

[Torsten]

Does anyone know why this is happening?

We must see your code to reproduce the problem.

[Phoebe Sayer]

@Torsten The equations are here

eqn1 = q1(t) + q2(t) - q7(t) == 0;
eqn2 = -q2(t) + q3(t) + q8(t) == 0;
eqn3 = -q1(t) -q3(t) + q9(t) == 0;
eqn4 = -q4(t) + q5(t) + q7(t) == 0;
eqn5 = -q5(t) + q6(t) - q8(t) == 0;
eqn6 = q4(t) -q6(t) - q9(t) == 0;
eqn7 = J*diff(q1(t)) + K*q1(t)^2 == p3(t) - p1(t);
eqn8 = J*diff(q2(t)) + K*q2(t)^2 == p2(t) - p1(t);
eqn9 = J*diff(q3(t)) + K*q3(t)^2 == p3(t) - p2(t);
eqn10 = J*diff(q4(t)) + K*q4(t)^2 == p4(t) - p6(t);
eqn11 = J*diff(q5(t)) + K*q5(t)^2 == p5(t) - p4(t);
eqn12 = J*diff(q6(t)) + K*q6(t)^2 == p6(t) - p5(t);
eqn13 = J*diff(q7(t)) + K*q7(t)^2 == p1(t) - p4(t);
eqn14 = J*diff(q8(t)) + K*q8(t)^2 == p5(t) - p2(t);
eqn15 = J*diff(q9(t)) + K*q9(t)^2 == p6(t) - p3(t);

Also of course p1 q1 etc were defined at the start with syms q1(t) q2(t) etc

And J and K are just numerical values that were defined at the start as well.

Answer 1

[Torsten]

The result of the following analysis of your differential-algebraic system is:

You must prescribe one of the six algebraic functions p1(t),p2(t),p3(t),p4(t),p5(t),p6(t) (e.g. p6(t) as shown below) because your algebraic relations between the q variables only have rank 5, not 6.

syms q1 q2 q3 q4 q5 q6 q7 q8 q9 
syms dq1 dq2 dq3 dq4 dq5 dq6 dq7 dq8 dq9 
syms p1 p2 p3 p4 p5 p6
J = 1;
K = 1;
eqn1 = q1 + q2 - q7 == 0;
eqn2 = -q2 + q3 + q8 == 0;
eqn3 = -q1 -q3 + q9 == 0;
eqn4 = -q4 + q5 + q7 == 0;
eqn5 = -q5 + q6 - q8 == 0;
eqn6 = q4 -q6 - q9 == 0;
[M,~] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6])

M = 

rank(M)

ans = 5

eqn1a = dq1 + dq2 - dq7 == 0;
eqn2a = -dq2 + dq3 + dq8 == 0;
eqn3a = -dq1 -dq3 + dq9 == 0;
eqn4a = -dq4 + dq5 + dq7 == 0;
eqn5a = -dq5 + dq6 - dq8 == 0;
eqn6a = dq4 -dq6 - dq9 == 0;
eqn7 = J*dq1 + K*q1^2 == p3 - p1;
eqn8 = J*dq2 + K*q2^2 == p2 - p1;
eqn9 = J*dq3 + K*q3^2 == p3 - p2;
eqn10 = J*dq4 + K*q4^2 == p4 - p6;
eqn11 = J*dq5 + K*q5^2 == p5 - p4;
eqn12 = J*dq6 + K*q6^2 == p6 - p5;
eqn13 = J*dq7 + K*q7^2 == p1 - p4;
eqn14 = J*dq8 + K*q8^2 == p5 - p2;
eqn15 = J*dq9 + K*q9^2 == p6 - p3;
s1 = solve([eqn1,eqn2,eqn4,eqn5,eqn6],[q1,q2,q4,q5,q9])

s1 = struct with fields:
    q1: q7 - q3 - q8
    q2: q3 + q8
    q4: q6 + q7 - q8
    q5: q6 - q8
    q9: q7 - q8

s1a = solve([eqn1a,eqn2a,eqn4a,eqn5a,eqn6a],[dq1,dq2,dq4,dq5,dq9])

s1a = struct with fields:
    dq1: dq7 - dq3 - dq8
    dq2: dq3 + dq8
    dq4: dq6 + dq7 - dq8
    dq5: dq6 - dq8
    dq9: dq7 - dq8

subs([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6],[q1,q2,q4,q5,q9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9])

ans = 

subs([eqn1a,eqn2a,eqn3a,eqn4a,eqn5a,eqn6a],[dq1,dq2,dq4,dq5,dq9],[s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

ans = 

eqn7 = subs(eqn7,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn7 = 

eqn8 = subs(eqn8,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn8 = 

eqn9 = subs(eqn9,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn9 = 

eqn10 = subs(eqn10,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn10 = 

eqn11 = subs(eqn11,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn11 = 

eqn12 = subs(eqn12,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn12 = 

eqn13 = subs(eqn13,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn13 = 

eqn14 = subs(eqn14,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn14 = 

eqn15 = subs(eqn15,[q1,q2,q4,q5,q9,dq1,dq2,dq4,dq5,dq9],[s1.q1,s1.q2,s1.q4,s1.q5,s1.q9,s1a.dq1,s1a.dq2,s1a.dq4,s1a.dq5,s1a.dq9])

eqn15 = 

solve([eqn7,eqn8,eqn9,eqn10,eqn11,eqn12,eqn13,eqn14,eqn15],[dq3,dq6,dq7,dq8,p1,p2,p3,p4,p5])

ans = struct with fields:
    dq3: (4*(q3 - q7 + q8)^2)/15 - (q7 - q8)^2/5 - (q6 - q8)^2/15 - (q6 + q7 - q8)^2/15 - (4*(q3 + q8)^2)/15 - (7*q3^2)/15 + (2*q6^2)/15 + q8^2/5
    dq6: (q7 - q8)^2/5 - (4*(q6 - q8)^2)/15 + (q3 - q7 + q8)^2/15 - (4*(q6 + q7 - q8)^2)/15 - (q3 + q8)^2/15 + (2*q3^2)/15 - (7*q6^2)/15 - q8^2/5
    dq7: (q6 - q8)^2/5 - (q7 - q8)^2/5 - (q3 - q7 + q8)^2/5 - (q6 + q7 - q8)^2/5 - (q3 + q8)^2/5 - (2*q7^2)/5 - q8^2/5
    dq8: (q6 - q8)^2/5 + (q7 - q8)^2/5 - (q3 + q8)^2/5 + q3^2/5 - q6^2/5 - q7^2/5 - (2*q8^2)/5
     p1: p6 - (q6 - q8)^2/15 - (2*(q7 - q8)^2)/5 - (q3 - q7 + q8)^2/3 + (q6 + q7 - q8)^2/3 - (4*(q3 + q8)^2)/15 - q3^2/15 - (4*q6^2)/15 + (2*q7^2)/5 - q8^2/5
     p2: p6 + (q6 - q8)^2/15 - (2*(q7 - q8)^2)/5 - (q3 - q7 + q8)^2/15 + (4*(q6 + q7 - q8)^2)/15 + (4*(q3 + q8)^2)/15 - q3^2/3 - q6^2/3 + q7^2/5 - (2*q8^2)/5
     p3: p6 - (3*(q7 - q8)^2)/5 + (q3 - q7 + q8)^2/5 + (q6 + q7 - q8)^2/5 + q3^2/5 - q6^2/5 + q7^2/5 - q8^2/5
     p4: p6 - (4*(q6 - q8)^2)/15 - (q7 - q8)^2/5 - (2*(q3 - q7 + q8)^2)/15 + (8*(q6 + q7 - q8)^2)/15 - (q3 + q8)^2/15 - q3^2/15 - (4*q6^2)/15 - q7^2/5
     p5: p6 + (4*(q6 - q8)^2)/15 - (q7 - q8)^2/5 - (q3 - q7 + q8)^2/15 + (4*(q6 + q7 - q8)^2)/15 + (q3 + q8)^2/15 - (2*q3^2)/15 - (8*q6^2)/15 + q8^2/5