solve the system of linear equation

I want to solve the system of linear equations, but one of these equations is a diffrential equation(the equation of Q), and one of these equations is described in if structure. I want to know the method of solving these, and my question is 'how can I solve this system'?

My simplified equations are written in below. My variables are Tg Tc Tb Ti Tw Q Tp and italic underlined alphabet are constant.

1)Tg=CTc

2)Tc=ETg+FTb

3)Tb=HTc+ITp

4)Ti=LTp

5)Tw=NTp

6)dQ/dt=OTb+PTp+QTi+RTw

7)if Q<Q_sol

Tp=DQ

elseif (Q_sol<Q) && (Q<Q_liq)

Tp=T_mp;

elseif Q > Q_liq

Tp=BQ;

end

[Walter Roberson]

That T_mp in the middle of the definition of T_p does not appear to be consistent with anything else??

Unfortunately the following approaches do not work (yet?)

syms Tb Tc Tg Ti Tp Tw T_mp
syms B_ubar C_ubar D_ubar E_ubar F_ubar H_ubar I_ubar L_ubar N_ubar O_ubar P_ubar Q_ubar R_ubar
syms Q(t) Q_sol Q_liq

eqn1 = Tg == C_ubar * Tc

eqn1 = 

eqn2 = Tc == E_ubar * Tg + F_ubar * Tb

eqn2 = 

eqn3 = Tb == H_ubar * Tc + I_ubar * Tp

eqn3 = 

eqn4 = Ti == L_ubar * Tp

eqn4 = 

eqn5 = Tw == N_ubar * Tp

eqn5 = 

dQ = diff(Q,t)

dQ(t) = 

eqn6 = dQ(t) == O_ubar * Tb + P_ubar * Tp + Q_ubar * Ti + R_ubar * Tw

eqn6 = 

Tp_defined = piecewise(Q(t) < Q_sol, D_ubar * Q(t), Q(t) < Q_liq, T_mp, B_ubar * Q(t))

Tp_defined = 

eqns = [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6].'

eqns = 

augmented = simplify(subs(eqns, Tp, Tp_defined))

augmented = 

sol = solve(augmented)

Warning: Unable to find explicit solution. For options, see help.

sol = struct with fields:
    T_mp: [0×1 sym]
      Tb: [0×1 sym]
      Tc: [0×1 sym]
      Tg: [0×1 sym]
      Ti: [0×1 sym]
       t: [0×1 sym]

Choose3 = @(var,LB,V1,UB,V2,V3) ((-V2 + V3)*heaviside(var - UB) - V1 + V2)*heaviside(var - LB) + V1

Choose3 = function_handle with value:
    @(var,LB,V1,UB,V2,V3)((-V2+V3)*heaviside(var-UB)-V1+V2)*heaviside(var-LB)+V1

Tp_defined2 = Choose3(Q(t), Q_sol, D_ubar * Q(t), Q_liq, T_mp, B_ubar * Q(t))

Tp_defined2 = 

augmented2 = simplify(subs(eqns, Tp, Tp_defined2))

augmented2 = 

sol2 = solve(augmented2)

Warning: Unable to find explicit solution. For options, see help.

sol2 = struct with fields:
    T_mp: [0×1 sym]
      Tb: [0×1 sym]
      Tc: [0×1 sym]
      Tg: [0×1 sym]
      Ti: [0×1 sym]
       t: [0×1 sym]

[Hannah Mohebi]

Thank you for your answer.T_mp is constant. I attached my real equations, and my code is below. My problem is with the choose3 function; I do not know how to write it for my equation; also, Matlab gives me a warning ' Unable to find an explicit solution.' How can I find the numerical solution?

clc;
clear;
close all;
syms T_b T_c T_g T_i T_p T_wmean Q_PCM(t)
I=50;
T_e=29.13+273;
vf=2;
vb=0;

%%%glass
alpha_g=0.05;             %Absorption coefficien of glass
delta_g=0.005;
landa_g=1;
h_ef=5.7+3.8*vf;
delta_c=0.0003;
landa_c=145;
h_ge=((delta_g/landa_g)+(1/h_ef))^(-1);
U_gc=((delta_g/landa_g)+(delta_c/landa_c))^(-1);
eqn1 = T_g ==(alpha_g*I+h_ge*T_e+U_gc*T_c)/(h_ge+U_gc);

%%%PVT
F=0.83;
alpha_c=0.9;
tha_g=0.91;
etha_ref=0.12;
beta=0.0045;
T_ref=293;
etha=(etha_ref*(1-(beta*(T_c-T_ref))));
delta_b=0.0001;
landa_b=0.36;
U_cb=((delta_c/landa_c)+(delta_b/landa_b))^(-1);
eqn2=T_c==(F*alpha_c*tha_g*I+U_gc*T_g-etha*F*alpha_c*tha_g*I+U_cb*T_b)/(U_gc+U_cb);

%%%Backplate
alpha_b=0.00039;
tha_c=0.09;
delta_p=0.05;
landa_PCMmean=0.195;
U_bp=((delta_b/landa_b)+(delta_p/landa_PCMmean))^(-1);
eqn3=T_b==(F*alpha_b*tha_c*tha_g*I+(1-F)*alpha_b*tha_g*I+U_cb*T_c+U_bp*T_p)/(U_bp+U_cb);

%%%insulatoin
delta_i=0.05;
landa_i=0.034;
h_eb=5.7+3.8*vb;
U_pi=((delta_i/landa_i)+(delta_p/landa_PCMmean))^(-1);
h_ie=((delta_i/landa_i)+(1/h_eb))^(-1);
eqn4 = T_i==(U_pi*T_p+h_ie*T_e)/(U_pi+h_ie);


%%%PCM layer
W=1;
L=1.8;
delta_t=0.0008;
landa_t=383;
delta_PCM=0.05;
landa_PCMl=0.15;
landa_PCMs=0.24;
Cp_PCMl=2100;
Cp_PCMs=2900;
rho_PCMl=780;
rho_PCMs=860;
LH=210000;
T_mp=28*273;
Q_sol=Cp_PCMs*(28-25);
Q_liq=Cp_PCMs*(28-25)+LH;
D=0.0254;
gama=1e-6;
c_w=4180;
mdot=0.03;
Dh=0.0254;
mu=1e-3;
landa_w=0.598;
A=pi*(D^2)/4;
Re=(mdot*Dh)/(mu*A);
Pr=c_w*mu/landa_w;
Nu=1.86*((Re*Pr*D/L)^(1/3));
h_pw=Nu*landa_w/D;

%%%water
T_win=298;
k_pw=((1/h_pw)+(delta_t/landa_t)+(delta_PCM/landa_PCMl))^(-1);
P=pi*D*h_pw/(c_w*mdot);
Q=exp(-P*L);
R=1-Q/L;
T_watL=T_win*Q+T_p*(1-Q);
eqn5 = T_wmean==T_win*R+T_p*(1-R); 

dQ_PCM = diff(Q_PCM,t);
eqn6 = dQ_PCM(t)==(W*L)*(U_bp*(T_b-T_p)-h_pw*(T_p-T_wmean)-U_pi*(T_p-T_i));
T_p_defined = piecewise(Q_PCM(t) < Q_sol, T_mp-(Q_sol-Q_PCM/Cp_PCMs), Q_PCM>Q_liq, T_mp-(Q_PCM-Q_liq/Cp_PCMl), T_mp);
eqns = [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6]';
augmented = simplify(subs(eqns, T_p, T_p_defined));
sol = solve(augmented);
Choose3 = @(var,LB,V1,UB,V2,V3) ((-V2 + V3)*heaviside(var - UB) - V1 + V2)*heaviside(var - LB) + V1;
T_p_defined2 = Choose3(Q_PCM(t), Q_sol, T_mp-(Q_sol-Q_PCM/Cp_PCMs), Q_liq, T_mp, T_mp-(Q_PCM-Q_liq/Cp_PCMl));
augmented2 = simplify(subs(eqns, T_p, T_p_defined2));
sol2 = solve(augmented2);

[Walter Roberson]

Choose3 = @(var,LB,V1,UB,V2,V3) is implementing

• var < LB --> output is V1
• var >= LB and var < UB --> output is V2
• var >= UB --> output is V3

in terms of heaviside instead of piecewise(), except that I do not promise that it is correct exactly at the boundaries (especially since I forget to set the HeavisideAtOrigin preference.)

Your image shows a three part test, and in each case it sets two output variables. The way to handle that in the context of using Choose3 is to use two different Choose3 calls, once for each output variable, with the var and LB and UB parameters the same for the two calls, but using different V1, V2, V3 values.

The reason I attempted with heaviside is that I know from experience that in some cases dsolve() can handle boundaries expressed as heaviside, but that it can rarely handle boundaries expressed as piecewise()

A problem is that we are invoking solve() rather than dsolve(). dsolve() fails on the linear parts, and solve() can fail on the differential parts.

It might be worth experimenting with breaking the system up into the three different conditions, getting out three systems that do not have internal conditions expressed in them. I would then follow that by seeing if it was possible to solve() the non-differential parts of that, hopefully leading to a single differential equation each that you would dsolve() .

When I look at your equations, it is not obvious to me that they are continuous equations -- not obvious to me that they have been chosen to match boundary conditions for the limit to be the same from both sides of the boundary.

Answer 1

[Torsten]

  Q0 = 356;
  tspan = [0,100];
  Qsol = 300;
  Qliq = 400;
  T_mp = ...;
  C = ...;
  E = ...;
  F = ...;
  H = ...;
  I = ...;
  L = ...;
  N = ...;
  O = ...;
  P = ...;
  Q = ...;
  R = ...;
  D = ...;
  B = ...;
  [t,Q] = ode45(@(t,Q)fun(t,Q,Qsol,Qliq,T_mp,C,E,F,H,I,L,N,O,P,Q,R,D,B),tspan,Q0);

function dQ = fun(t,Q,Qsol,Qliq,T_mp,C,E,F,H,I,L,N,O,P,Q,R,D,B) 
  if Q < Qsol
    Tp = D*Q;
  elseif Q >= Qsol && Q <= Qliq
    Tp = T_mp;
  elseif Q > Qliq
    Tp = B*Q;
  end
  Tw = N*Tp;
  Ti = L*Tp;
  Tc = F*I/(-F*H+1-E*C)*Tp;
  Tg = C*Tc;
  Tb = I*Tp+H*Tc;
  dQ = O*Tb+P*Tp+Q*Ti+R*Tw;
end

[Walter Roberson]

  if Q < Qsol
    Tp = D*Q;
  elseif Q >= Qsol && Q <= Qliq
    Tp = T_mp;
  elseif Q > Qliq
    Tp = B*Q;
  end

That is likely to either fail or give wrong answers. Using if inside of an ode function only works mathematically under one of two conditions:

1. Only a single branch of the if is used for any one ode*() function call; or
2. The branches have direct values at the boundaries that are equal, and first derivatives at the boundaries that are equal, and second derivatives at the boundaries that are equal. continuity at least.

Because the value being tested, Q, is not time-based, it would probably be necessary to use event functions to detect crossing a boundary, and have the system stop integration at that point; you would then call ode45() again to continue on the other side of the boundary. Sounds pointless, I know, but the ode*() functions are doing internal interpolation, and that interpolation becomes invalid each time you cross a boundary.