MATLAB Answers

Dealing with symbolic parameter in a DAE

5 views (last 30 days)
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.

  2 Comments

Christopher Lamb
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);

Sign in to comment.

Accepted Answer

Guru Mohanty
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
ode15s_case.jpg

  1 Comment

Sign in to comment.

More Answers (0)

Sign in to answer this question.


Translated by