error [empty sym] using dsolve "Warning: unable to find symbolic solution"
8 views (last 30 days)
Show older comments
hi. i'm trying to create a simutation program for a refrigerator and i'm having problems.
i'm trying to solve the differntioal equations below but MATLAB returns me an error [empty sym] "Warning: unable to find symbolic solution".
if i have anything wrong in the equations or if someone knows why i'm having this error, it would be very thankful if you let me know.
and please excuse my poor english (i'm an student from japan).
syms M_w c_pco c_pev c_phx c_pr c_ps c_pw D D_s k_a L_w m_c_co m_c_hx m_ch M_co M_ev m_h M_hx M_s p_r q_max Q_s r_s T_c_co_in T_c_hx_in T_ch_in T_h_in h_r
syms t q_de(t) T_de(t) q_ad(t) T_ad(t) M_w_co(t) T_co(t) T_ev(t) M_w_ev(t) T_ch_out(t) T_h_out(t) T_c_hx_out(t) T_c_co_out(t) q_0_ad(t) q_0_de(t)
M_s = 0.69
M_w = 0.789
c_ps = 0.84
c_pr = 2.5*10^(-3)
c_pw = 0.42
M_hx = 4.6282
c_phx = 0.5
m_h = 0.333
T_h_in = 353
m_c_hx = 0.267
T_c_hx_in = 303
M_co = 6.0547
c_pco = 0.50
L_w = 2.42
m_c_co = 0.267
T_c_co_in = 303
M_ev = 6.0547
c_pev = 0.50
m_ch = 0.1
T_ch_in = 300
D_s = 0.63*10^(-10)
k_a = 1.0
q_max = 1.24
D = 2.844*10^(-6)
p_r = 0.56
r_s = 2*10^(-3)
Q_s = 1.2715
h_r = 38.6/46.07*1000
eqn1(t) = (M_s*c_ps+M_s*c_pr*q_de+M_hx*c_phx)*diff(T_de, t) == m_h*c_pw*(T_h_in-T_h_out)+Q_s*diff(q_de, t)
eqn2(t) = (M_s*c_ps+M_s*c_pr*q_ad+M_hx*c_phx)*diff(T_ad, t) == m_h*c_pw*(T_c_hx_in-T_c_hx_out)+Q_s*diff(q_ad, t)
eqn3(t) = (M_w_co*c_pw+M_co*c_pco)*diff(T_co, t) == -h_r*diff(q_de, t)+c_pr*(T_de-T_co)*diff(q_de, t)+m_c_co*c_pw*(T_c_co_in-T_c_co_out)
eqn4(t) = (M_w_ev*c_pw+M_ev*c_pev)*diff(T_ev, t) == -h_r*diff(q_ad, t)+m_ch*c_pw*(T_ch_in-T_ch_out)
eqn5_de(t) = T_h_out == T_de+(T_h_in-T_de)*exp(-0.88)
eqn5_ad(t) = T_c_hx_out == T_ad+(T_c_hx_in-T_ad)*exp(-0.88)
eqn5_co(t) = T_c_co_out == T_co+(T_c_co_in-T_co)*exp(-0.88)
eqn5_ev(t) = T_ch_out == T_ev+(T_ch_in-T_ev)*exp(-0.88)
eqn6_1(t) = diff(M_w_co, t)+diff(q_de, t) == 0
eqn6_2(t) = diff(M_w_ev, t)+diff(q_ad, t) == 0
eqn7_ad(t) = q_0_ad == q_max*exp(-D*(T_ad*log(p_r))^2)
eqn8_ad(t) = diff(q_ad, t) == 15*D_s\r_s^2*(q_0_ad-q_ad)
eqn7_de(t) = q_0_de == q_max*exp(-D*(T_de*log(p_r))^2)
eqn8_de(t) = diff(q_de, t) == 15*D_s\r_s^2*(q_0_de-q_de)
eqn1 = subs(eqn1, lhs(eqn5_de), rhs(eqn5_de))
eqn2 = subs(eqn2, lhs(eqn5_ad), rhs(eqn5_ad))
eqn3 = subs(eqn3, lhs(eqn5_co), rhs(eqn5_co))
eqn4 = subs(eqn4, lhs(eqn5_ev), rhs(eqn5_ev))
eqn8_ad = subs(eqn8_ad, lhs(eqn7_ad), rhs(eqn7_ad))
eqn8_de = subs(eqn8_de, lhs(eqn7_de), rhs(eqn7_de))
eqns = [eqn1(t); eqn2(t); eqn3(t); eqn4(t); eqn6_1(t); eqn6_2(t); eqn8_ad(t); eqn8_de(t)]
cond_T_ad = T_ad(0) == T_c_hx_in
cond_T_de = T_de(0) == T_h_in
cond_T_co = T_co(0) == T_c_co_in
cond_T_ev = T_ev(0) == T_ch_in
cond_T_c_co_out = T_c_co_out(0) == T_c_co_in
cond_T_c_hx_out = T_c_hx_out(0) == T_c_hx_in
cond_T_ch_out = T_ch_out(0) == T_ch_in
cond_T_h_out = T_h_out(0) == T_h_in
cond_q_ad = q_ad(0) == 0
cond_q_de = q_de(0) == 0
cond_M_w_co = M_w_co(0) == 0
cond_M_w_ev = M_w_ev(0) == M_w
conds = [cond_T_ad; cond_T_de; cond_T_co; cond_T_ev; cond_q_ad; cond_q_de; cond_M_w_co; cond_M_w_ev]
S = dsolve(eqns, conds)
0 Comments
Answers (3)
Star Strider
on 26 Nov 2021
It is likely not possible to derive a closed-form (analytic or symbolic) solution for those. Create an anonymous function to use with one of the numerical oODE solvers instead, and integrate with one of them.
syms M_w c_pco c_pev c_phx c_pr c_ps c_pw D D_s k_a L_w m_c_co m_c_hx m_ch M_co M_ev m_h M_hx M_s p_r q_max Q_s r_s T_c_co_in T_c_hx_in T_ch_in T_h_in h_r
syms t q_de(t) T_de(t) q_ad(t) T_ad(t) M_w_co(t) T_co(t) T_ev(t) M_w_ev(t) T_ch_out(t) T_h_out(t) T_c_hx_out(t) T_c_co_out(t) q_0_ad(t) q_0_de(t) Y
sympref('AbbreviateOutput',false);
M_s = 0.69
M_w = 0.789
c_ps = 0.84
c_pr = 2.5*10^(-3)
c_pw = 0.42
M_hx = 4.6282
c_phx = 0.5
m_h = 0.333
T_h_in = 353
m_c_hx = 0.267
T_c_hx_in = 303
M_co = 6.0547
c_pco = 0.50
L_w = 2.42
m_c_co = 0.267
T_c_co_in = 303
M_ev = 6.0547
c_pev = 0.50
m_ch = 0.1
T_ch_in = 300
D_s = 0.63*10^(-10)
k_a = 1.0
q_max = 1.24
D = 2.844*10^(-6)
p_r = 0.56
r_s = 2*10^(-3)
Q_s = 1.2715
h_r = 38.6/46.07*1000
eqn1(t) = (M_s*c_ps+M_s*c_pr*q_de+M_hx*c_phx)*diff(T_de, t) == m_h*c_pw*(T_h_in-T_h_out)+Q_s*diff(q_de, t)
eqn2(t) = (M_s*c_ps+M_s*c_pr*q_ad+M_hx*c_phx)*diff(T_ad, t) == m_h*c_pw*(T_c_hx_in-T_c_hx_out)+Q_s*diff(q_ad, t)
eqn3(t) = (M_w_co*c_pw+M_co*c_pco)*diff(T_co, t) == -h_r*diff(q_de, t)+c_pr*(T_de-T_co)*diff(q_de, t)+m_c_co*c_pw*(T_c_co_in-T_c_co_out)
eqn4(t) = (M_w_ev*c_pw+M_ev*c_pev)*diff(T_ev, t) == -h_r*diff(q_ad, t)+m_ch*c_pw*(T_ch_in-T_ch_out)
eqn5_de(t) = T_h_out == T_de+(T_h_in-T_de)*exp(-0.88)
eqn5_ad(t) = T_c_hx_out == T_ad+(T_c_hx_in-T_ad)*exp(-0.88)
eqn5_co(t) = T_c_co_out == T_co+(T_c_co_in-T_co)*exp(-0.88)
eqn5_ev(t) = T_ch_out == T_ev+(T_ch_in-T_ev)*exp(-0.88)
eqn6_1(t) = diff(M_w_co, t)+diff(q_de, t) == 0
eqn6_2(t) = diff(M_w_ev, t)+diff(q_ad, t) == 0
eqn7_ad(t) = q_0_ad == q_max*exp(-D*(T_ad*log(p_r))^2)
eqn8_ad(t) = diff(q_ad, t) == 15*D_s\r_s^2*(q_0_ad-q_ad)
eqn7_de(t) = q_0_de == q_max*exp(-D*(T_de*log(p_r))^2)
eqn8_de(t) = diff(q_de, t) == 15*D_s\r_s^2*(q_0_de-q_de)
eqn1 = subs(eqn1, lhs(eqn5_de), rhs(eqn5_de))
eqn2 = subs(eqn2, lhs(eqn5_ad), rhs(eqn5_ad))
eqn3 = subs(eqn3, lhs(eqn5_co), rhs(eqn5_co))
eqn4 = subs(eqn4, lhs(eqn5_ev), rhs(eqn5_ev))
eqn8_ad = subs(eqn8_ad, lhs(eqn7_ad), rhs(eqn7_ad))
eqn8_de = subs(eqn8_de, lhs(eqn7_de), rhs(eqn7_de))
eqns = [eqn1(t); eqn2(t); eqn3(t); eqn4(t); eqn6_1(t); eqn6_2(t); eqn8_ad(t); eqn8_de(t)]
% cond_T_ad = T_ad(0) == T_c_hx_in
% cond_T_de = T_de(0) == T_h_in
% cond_T_co = T_co(0) == T_c_co_in
% cond_T_ev = T_ev(0) == T_ch_in
% cond_T_c_co_out = T_c_co_out(0) == T_c_co_in
% cond_T_c_hx_out = T_c_hx_out(0) == T_c_hx_in
% cond_T_ch_out = T_ch_out(0) == T_ch_in
% cond_T_h_out = T_h_out(0) == T_h_in
% cond_q_ad = q_ad(0) == 0
% cond_q_de = q_de(0) == 0
% cond_M_w_co = M_w_co(0) == 0
% cond_M_w_ev = M_w_ev(0) == M_w
% conds = [cond_T_ad; cond_T_de; cond_T_co; cond_T_ev; cond_q_ad; cond_q_de; cond_M_w_co; cond_M_w_ev]
% S = dsolve(eqns, conds)
[VF,Sbs] = odeToVectorField(eqns)
eqnsfcn = matlabFunction(VF, 'Vars',{'t','Y','q_de','T_de','T_ad','q_ad','T_co','T_ev','M_w_co','M_w_ev'})
Then use ode45 (or other suitable solver, perhaps ode15s if the system turns out to be stiff. (If ode45 takes longer than a minute or two to integrate these, ithe system is likely stiff. Just change the solver to ode15s or one of the other stiff solvers, and it should integrate quickly.)
The ode45 call should be something like this —
initcons = [ ... ]; % Initial Conditions Vector
start_time = 0;
end_time = ...;
tspan = [start_time end_time];
[t,y] = ode45(@(t,y)eqnsfcn(t,Y,q_de,T_de,T_ad,q_ad,T_co,T_ev,M_w_co,M_w_ev), tspan, initconds);
figure
plot(t, y)
grid
.
12 Comments
Walter Roberson
on 29 Dec 2021
Your revised system of equations does not take the derivative of T_ad(t), T_co(t), T_de(t), or T_ev(t) .
Walter Roberson
on 30 Nov 2021
conds = [cond_T_ad; cond_T_de; cond_T_co; cond_T_ev; cond_q_ad; cond_q_de; cond_M_w_co; cond_M_w_ev]
Those are all equations
initconds = [conds]
and you set initconds to that
[t, y] = ode45(@(t, y)eqnsfcn(t, Y, q_de, T_de, T_ad, q_ad, T_co, T_ev, M_w_co, M_w_ev), tspan, initconds)
and you pass the equations to ode45. But ode45 needs numeric initial conditions.
At any point, the boundary conditions for the ode are going to be passed into the function as the second parameter, where they are received as y (lower-case) . But you never pass y (lower-case) into your function: you pass Y (upper-case), which is a scalar symbolic variable.
0 Comments
See Also
Categories
Find more on Symbolic Math Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!