Dealing with symbolic parameter in a DAE

Hello!

I've been using the Mass Matrix Solvers page to create a ODE solver for a mass matrix that I'm using for a piston.

However, I've come across a hickup where I get an error due to a symbolic parameter being in the DAE eqn.

Here since a screenshot of the matrix:

On the left matrix you can see in the third row, the specific heat is the parameter that I mentioned.

Here is the error I'm getting:

Error using mupadengine/feval (line 187)
Found symbolic object 'Cv' in DAEs. Only variables and declared parameters can be symbolic

What I'm asking is how am I able to create a work around for this.

[Star Strider]

Code?

[Christopher Lamb]

syms p(theta) m(theta) T(theta) W(theta) Q(theta) V(theta) Cv mdotin...
      mdotex mdotleak Cp Tair Qdotrxn Qdotloss omega B L a

eqn1 = diff(p(theta), 1)/p(theta)-diff(m(theta), 1)/m(theta)...
        -diff(T(theta),1)/T(theta)+diff(V(theta), 1)/V(theta) == 0;
eqn2 = diff(m(theta), 1) == 1/omega*(mdotin - mdotex - mdotleak);
eqn3 = Cv*T(theta)*diff(m(theta), 1)+Cv*m(theta)*diff(T(theta), 1)...
        -diff(Q(theta), 1)+diff(W(theta), 1) == 0;
eqn4 = diff(W(theta), 1)-p(theta)*diff(V(theta), 1) == 0;
eqn5 = diff(Q(theta), 1) == 1/omega*(mdotin*Cp*Tair-mdotex*Cp*T(theta)...
        -mdotleak*Cp*T(theta)+Qdotrxn-Qdotloss);
eqn6 = diff(V(theta), 1) == 1/omega*((B^2*pi*(a*sin(theta) +...
        (a^2*cos(theta)*sin(theta))/(- a^2*sin(theta)^2 + L^2)^(1/2)))/4);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6];

vars = [p(theta); m(theta); T(theta); W(theta); Q(theta); V(theta)];
origVars = length(vars);

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)

[M,f] = massMatrixForm(DAEs,DAEvars)

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars)

M = odeFunction(M, DAEvars);
f = odeFunction(f, DAEvars,Cv,mdotin,mdotex,mdotleak,Cp,Tair,...
                 Qdotrxn,Qdotloss,omega, B, L, a);


Cv = 0.718;
mdotin = 0;
mdotex = 0;
mdotleak = 0;
Cp = 1.005;
Tair = 273+21;
Qdotrxn = 0;
Qdotloss = 0;
B = 86;
L = 50;
a = 150;
F = @(theta, Y) f(theta, Y, Cv, mdotin, mdotex, mdotleak, Cp, Tair, Qdotrxn,...
              Qdotloss, omega, B, L, a);

Answer 1

[Guru Mohanty]

Hi, I understand you are trying to solve DAEs Using Mass Matrix Solvers. The error is due to missing input argument of the odeFunction. However, I can get solutions to these DAEs considering zero initial condition. Here is the code for it.

clc;
clear all;
syms p(theta) m(theta) T(theta) W(theta) Q(theta) V(theta) Cv mdotin...
      mdotex mdotleak Cp Tair Qdotrxn Qdotloss  B L a
eqn1 = diff(p(theta), 1)/p(theta)-diff(m(theta), 1)/m(theta)...
        -diff(T(theta),1)/T(theta)+diff(V(theta), 1)/V(theta) == 0;
eqn2 = diff(m(theta), 1) == (mdotin - mdotex - mdotleak);
eqn3 = Cv*T(theta)*diff(m(theta), 1)+Cv*m(theta)*diff(T(theta), 1)...
        -diff(Q(theta), 1)+diff(W(theta), 1) == 0;
eqn4 = diff(W(theta), 1)-p(theta)*diff(V(theta), 1) == 0;
eqn5 = diff(Q(theta), 1) == (mdotin*Cp*Tair-mdotex*Cp*T(theta)...
        -mdotleak*Cp*T(theta)+Qdotrxn-Qdotloss);
eqn6 = diff(V(theta), 1) == ((B^2*pi*(a*sin(theta) +...
        (a^2*cos(theta)*sin(theta))/(- a^2*sin(theta)^2 + L^2)^(1/2)))/4);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6];
vars = [p(theta); m(theta); T(theta); W(theta); Q(theta); V(theta)];
origVars = length(vars);
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
[M,f] = massMatrixForm(DAEs,DAEvars);
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars);
M = odeFunction(M, DAEvars,Cv);
f = odeFunction(f, DAEvars,Cv,mdotin,mdotex,mdotleak,Cp,Tair,...
                 Qdotrxn,Qdotloss, B, L, a);
Cv = 0.718;
mdotin = 0;
mdotex = 0;
mdotleak = 0;
Cp = 1.005;
Tair = 273+21;
Qdotrxn = 0;
Qdotloss = 0;
B = 86;
L = 50;
a = 150;
F = @(theta, Y) f(theta, Y, Cv, mdotin, mdotex, mdotleak, Cp, Tair, Qdotrxn,...
              Qdotloss, B, L, a);
% Zero Initial Condition
y0=zeros(6,1);
% Solve System of ODE
[tSol,ySol] = ode15s(F, [0, 0.5], y0,0);
plot(tSol,ySol(:,1:origVars),'-o')
 
for k = 1:origVars
  S{k} = char(DAEvars(k));
end
grid on

[Christopher Lamb]

Thank you for the help!