# Dealing with symbolic parameter in a DAE

5 views (last 30 days)
Christopher Lamb on 15 Jan 2020
Commented: Christopher Lamb on 29 Jan 2020
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:
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 on 15 Jan 2020
Code?
Christopher Lamb on 16 Jan 2020
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);

Guru Mohanty on 23 Jan 2020
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

#### 1 Comment

Christopher Lamb on 29 Jan 2020
Thank you for the help!