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

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)

Answer 1

[Star Strider]

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

.

[Ko Nakahara]

thank you very much.

I entered the initconds and end time and ran the program, but it says "index exceeds the number of array elements (1)" when I used the ode45.

I'm sorry to take your time, but it would be very thankful if you let me know what this error means.

[Star Strider]

My pleasure.

Since I do not have the rest of the code being used, specifically what ‘initcons’ is, I cannot run the code myself to check it, so I have no idea what could be throwing that error. There are 8 differential equations, so there should be 8 initial conditions specified.

Since matlabFunction created the anonymous function, sometimes this is caused by using a row vector when an column vector is appropriate, or the other way round. I have no way of knowing without seeing the rest of the code being used for the ode45 call.

.

[Star Strider]

There were two problems.

First, ‘initconds’ must be defined as:

initconds = double(rhs(conds));

and there was a typographical error in the ode45 call.

It should be:

[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);

With those two corrections, it runs without error (although the output is not very exciting).

The full calling sequence (after all the symbolic results) is then:

[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 = double(rhs(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
loglog(t, y)
grid
legend(string(sbs), 'Location','best')

So with those changes, it works, as the plot demonstrates.

.

[Ko Nakahara]

thank you very much.

I managed to run the code, but I have one more question.

there are 8 data in the graph, but is it possible to plot only specific data? (for example, just plot q_ad and q_de and compare them)

sorry for asking basic questions.

[Walter Roberson]

Select columns from y

plot(t, y(:, [3, 6, 7]))

changing the numbers to the column number you need

[Star Strider]

@Ko Nakahara — My pleasure!

No worries about the questions.

I was sleeping when this Comment was posted, and Walter responded to it. I am not certain what ‘compare’ means in this context, however they were all calculated at the same time instants, so that sho8uld not be difficult.

Note that there are only about four distinct lines on the plot, so since I did not restrict them in any way, some variables may be over-plotting others. (I did not check specifically to see which ones they are.)

@Walter Roberson — Thank you.

.

[Ko Nakahara]

@Star Strider @Walter Roberson

thank you for the help.

I realized there were a few errors in the equation, so I fixed it.

("diff(T_x, t)" in the left side of each equations from 1 through 4 was actually "T_x")

then when I used "[VF, sbs] = odeToVectorField(eqns)" in line 66, it didn't work.

does this mean I have to use a different method?

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)*T_de == 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)*T_ad == 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)*T_co == -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)*T_ev == -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_h_in
cond_T_de = T_de(0) == T_c_hx_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_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)
figure
plot(t, y(:, [1, 4]))
grid

[Star Strider]

‘does this mean I have to use a different method?’

Yes.

[Ko Nakahara]

@Star Strider

thank you for the response.

i have been working on it, but i haven't been able to run the program.

could you give me some advise on how to run it?

[Star Strider]

I honestly have no idea what the code does.

Create the differential equations as an anonymous function that the numeric differential equation solvers can run, as I previously demonstrated.

The functions (that were previously part of the differential equation system) need to be transformed into anonymous functions themselves, so the differential equations can call them as necessary.

The initial conditions must be created as a vector that the numerical differential equation solvers can use.

That is the only advise I can provide.

[Walter Roberson]

Your revised system of equations does not take the derivative of T_ad(t), T_co(t), T_de(t), or T_ev(t) .

[Ko Nakahara]

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)

Answer 2

[Ko Nakahara]

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

[Nikhil Kapoor]

How to get values of q_de, T_de, T_ad etc

Answer 3

[Walter Roberson]

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.