error [empty sym] using dsolve "Warning: unable to find symbolic solution"

19 views (last 30 days)
Ko Nakahara
Ko Nakahara on 26 Nov 2021
Commented: Ko Nakahara on 26 Jan 2022 at 12:48
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)

Answers (3)

Star Strider
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_s = 0.6900
M_w = 0.789
M_w = 0.7890
c_ps = 0.84
c_ps = 0.8400
c_pr = 2.5*10^(-3)
c_pr = 0.0025
c_pw = 0.42
c_pw = 0.4200
M_hx = 4.6282
M_hx = 4.6282
c_phx = 0.5
c_phx = 0.5000
m_h = 0.333
m_h = 0.3330
T_h_in = 353
T_h_in = 353
m_c_hx = 0.267
m_c_hx = 0.2670
T_c_hx_in = 303
T_c_hx_in = 303
M_co = 6.0547
M_co = 6.0547
c_pco = 0.50
c_pco = 0.5000
L_w = 2.42
L_w = 2.4200
m_c_co = 0.267
m_c_co = 0.2670
T_c_co_in = 303
T_c_co_in = 303
M_ev = 6.0547
M_ev = 6.0547
c_pev = 0.50
c_pev = 0.5000
m_ch = 0.1
m_ch = 0.1000
T_ch_in = 300
T_ch_in = 300
D_s = 0.63*10^(-10)
D_s = 6.3000e-11
k_a = 1.0
k_a = 1
q_max = 1.24
q_max = 1.2400
D = 2.844*10^(-6)
D = 2.8440e-06
p_r = 0.56
p_r = 0.5600
r_s = 2*10^(-3)
r_s = 0.0020
Q_s = 1.2715
Q_s = 1.2715
h_r = 38.6/46.07*1000
h_r = 837.8554
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)
eqn1(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)
eqn2(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)
eqn3(t) = 
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)
eqn4(t) = 
eqn5_de(t) = T_h_out == T_de+(T_h_in-T_de)*exp(-0.88)
eqn5_de(t) = 
eqn5_ad(t) = T_c_hx_out == T_ad+(T_c_hx_in-T_ad)*exp(-0.88)
eqn5_ad(t) = 
eqn5_co(t) = T_c_co_out == T_co+(T_c_co_in-T_co)*exp(-0.88)
eqn5_co(t) = 
eqn5_ev(t) = T_ch_out == T_ev+(T_ch_in-T_ev)*exp(-0.88)
eqn5_ev(t) = 
eqn6_1(t) = diff(M_w_co, t)+diff(q_de, t) == 0
eqn6_1(t) = 
eqn6_2(t) = diff(M_w_ev, t)+diff(q_ad, t) == 0
eqn6_2(t) = 
eqn7_ad(t) = q_0_ad == q_max*exp(-D*(T_ad*log(p_r))^2)
eqn7_ad(t) = 
eqn8_ad(t) = diff(q_ad, t) == 15*D_s\r_s^2*(q_0_ad-q_ad)
eqn8_ad(t) = 
eqn7_de(t) = q_0_de == q_max*exp(-D*(T_de*log(p_r))^2)
eqn7_de(t) = 
eqn8_de(t) = diff(q_de, t) == 15*D_s\r_s^2*(q_0_de-q_de)
eqn8_de(t) = 
eqn1 = subs(eqn1, lhs(eqn5_de), rhs(eqn5_de))
eqn1(t) = 
eqn2 = subs(eqn2, lhs(eqn5_ad), rhs(eqn5_ad))
eqn2(t) = 
eqn3 = subs(eqn3, lhs(eqn5_co), rhs(eqn5_co))
eqn3(t) = 
eqn4 = subs(eqn4, lhs(eqn5_ev), rhs(eqn5_ev))
eqn4(t) = 
eqn8_ad = subs(eqn8_ad, lhs(eqn7_ad), rhs(eqn7_ad))
eqn8_ad(t) = 
eqn8_de = subs(eqn8_de, lhs(eqn7_de), rhs(eqn7_de))
eqn8_de(t) = 
eqns = [eqn1(t); eqn2(t); eqn3(t); eqn4(t); eqn6_1(t); eqn6_2(t); eqn8_ad(t); eqn8_de(t)]
eqns = 
% 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)
VF = 
Sbs = 
eqnsfcn = matlabFunction(VF, 'Vars',{'t','Y','q_de','T_de','T_ad','q_ad','T_co','T_ev','M_w_co','M_w_ev'})
eqnsfcn = function_handle with value:
@(t,Y,q_de,T_de,T_ad,q_ad,T_co,T_ev,M_w_co,M_w_ev)[exp(Y(2).^2.*(-9.561229022719256e-7)).*(exp(Y(2).^2.*9.561229022719256e-7).*Y(1).*2.5e+1-3.1e+1).*(-1.693121693121693e+2);(exp(Y(2).^2.*(-9.561229022719256e-7)).*(exp(Y(2).^2.*9.561229022719256e-7).*1.643970740393376e+19-exp(Y(2).^2.*9.561229022719256e-7).*Y(1).*3.062339968622479e+21-exp(Y(2).^2.*9.561229022719256e-7).*Y(2).*4.65714090762996e+16+3.797301561091874e+21).*7.029932191926491e-14)./(Y(1).*6.9e+1+1.15748e+5);(exp(Y(3).^2.*(-9.561229022719256e-7)).*(exp(Y(3).^2.*9.561229022719256e-7).*1.411113695011878e+19-exp(Y(3).^2.*9.561229022719256e-7).*Y(3).*4.65714090762996e+16-exp(Y(3).^2.*9.561229022719256e-7).*Y(4).*3.062339968622479e+21+3.797301561091874e+21).*7.029932191926491e-14)./(Y(4).*6.9e+1+1.15748e+5);exp(Y(3).^2.*(-9.561229022719256e-7)).*(exp(Y(3).^2.*9.561229022719256e-7).*Y(4).*2.5e+1-3.1e+1).*(-1.693121693121693e+2);(exp(Y(2).^2.*(-9.561229022719256e-7)).*(exp(Y(2).^2.*9.561229022719256e-7).*3.796460854590528e+26+exp(Y(2).^2.*9.561229022719256e-7).*Y(1).*6.771050021228565e+31-exp(Y(2).^2.*9.561229022719256e-7).*Y(5).*1.25295737775265e+24+Y(2).*2.505235883113471e+26-Y(5).*2.505235883113471e+26-exp(Y(2).^2.*9.561229022719256e-7).*Y(1).*Y(2).*2.020351518639896e+26+exp(Y(2).^2.*9.561229022719256e-7).*Y(1).*Y(5).*2.020351518639896e+26-8.396102026323421e+31).*1.047541527737154e-21)./(Y(7).*8.4e+3+6.0547e+4);(exp(Y(3).^2.*(-9.561229022719256e-7)).*(exp(Y(3).^2.*9.561229022719256e-7).*8.044670162135792e+23+exp(Y(3).^2.*9.561229022719256e-7).*Y(4).*3.869171440702037e+29-exp(Y(3).^2.*9.561229022719256e-7).*Y(6).*2.681556720711931e+21-4.797772586470526e+29).*1.83319767354002e-19)./(Y(8).*8.4e+3+6.0547e+4);exp(Y(2).^2.*(-9.561229022719256e-7)).*(exp(Y(2).^2.*9.561229022719256e-7).*Y(1).*2.5e+1-3.1e+1).*1.693121693121693e+2;exp(Y(3).^2.*(-9.561229022719256e-7)).*(exp(Y(3).^2.*9.561229022719256e-7).*Y(4).*2.5e+1-3.1e+1).*1.693121693121693e+2]
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
Ko Nakahara
Ko Nakahara on 26 Jan 2022 at 12:48
sorry for the late response.
i have managed to run the program i was working on by the script below.
thank you for your help.
now i am trying to output the plot of Q_ch in eqn9_ch.
can you give me some advise on how to do it?
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_ad(t) q_de(t) T_ad(t) T_de(t) T_co(t) T_ev(t) M_w_co(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 Q_ch(t) Q_ad(t) Q_h(t)
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 = 17/60
T_c_hx_in = 303
M_co = 6.0547
c_pco = 0.50
L_w = 2.42
m_c_co = 17/60
T_c_co_in = 303
M_ev = 6.0547
c_pev = 0.50
m_ch = 17/60
T_ch_in = 300
D_s = 0.63*10^(-10)
k_a = 1.0
q_max = 0.5
D = 2.844*10^(-6)
p_r_ad = 1/0.56
p_r_de= 1/0.11
r_s = 2*10^(-3)
Q_s = 29.9246
h_r = 38.6/46.07*1000
eqn1(t) = (M_s*c_ps+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+c_pr*q_ad+M_hx*c_phx)*diff(T_ad, t) == m_c_hx*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)+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_ad))^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_de))^2)
eqn8_de(t) = diff(q_de, t) == 15*D_s\r_s^2*(q_0_de-q_de)
eqn9_ch(t) = Q_ch(t) == m_ch*c_pw*(T_ch_in-T_ch_out)
eqn9_ad(t) = Q_ad(t) == m_c_hx*c_pw*(T_c_hx_out-T_c_hx_in)
eqn9_h(t) = Q_h(t) == m_h*c_pw*(T_h_in-T_h_out)
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_q_de = q_de(0) == 0.5
cond_T_de = T_de(0) == T_h_in
cond_T_ad = T_ad(0) == T_c_hx_in
cond_q_ad = q_ad(0) == 0
cond_T_co = T_co(0) == T_c_co_in
cond_T_ev = T_ev(0) == T_ch_in
cond_M_w_co = M_w_co(0) == 0
cond_M_w_ev = M_w_ev(0) == 0.789
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
conds = [cond_q_de; cond_T_de; cond_T_ad; cond_q_ad; cond_T_co; cond_T_ev; cond_M_w_co; cond_M_w_ev]
[VF, sbs] = odeToVectorField(eqns)
eqnsfcn = matlabFunction(VF, 'Vars', {'t', 'Y', 'q_ad', 'q_de', 'T_ad', 'T_de', 'T_co', 'T_ev', 'M_w_co', 'M_w_ev'})
initconds = double(rhs(conds))
t_start = 0
t_end = 400
tspan = [t_start t_end]
[t, y] = ode45(@(t, y)eqnsfcn(t, y, q_ad, q_de, T_ad, T_de, T_co, T_ev, M_w_co, M_w_ev), tspan, initconds)

Sign in to comment.


Ko Nakahara
Ko Nakahara on 30 Nov 2021
Sorry for the lack of explanation.
I used the 8 initial conditions in "conds" for "initconds".
I attached the full code below.
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]
[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'})
initconds = [conds]
t_start = 0
t_end = 400
tspan = [t_start t_end]
[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

Walter Roberson
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.

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!