Dear Maple enthusiasts,

I am unable to find a working method to solve a system of 8 equations, of which 4 are differential equations. The system contains 8 unknown variables and the goal is to find an expression for each of these variables as a function of the time t. I have attached the code of my project at the bottom of this message.

I have tried the following:

1. Using solve/dsolve to solve all 8 equations at once. This results in Maple eating up all of my memory and never finishing its calculations.
2. First using solve to solve the 4 non-differential equations so that I get 4 out of 8 variables as a function of the 4 remaining variables. This results in an expression containing RootOf() for each of the 4 veriables I'm solving for, which prevents me from using these expressions in the 4 remaining differential equations.
3. First using dsolve to solve the differential equations, which gives once again an expression for 4 variables as a function of the 4 remaining variables. I then use solve to solve the 4 remaining equations with the new found expressions. This results in an extremely long solution for each of the variables.

The code below contains the 3rd option I tried.

Any help or suggestions would be greatly appreciated. I have been scratching my head so much that I'm getting bald and whatever I search for on google or in the Maple help, I can't find a good reference to a system of differential equations together with other equations.

Download Park_-_mixed_control_model.mw

Dear Maple enthusiasts,

I was debugging a piece of Maple code I wrote with 3 for loops in it, because Maple kept running infinitely when I executed the worksheet. I found out that in the inner most loop Maple was getting stuck on integrating a function with a discontinuity. However, at the same time I noticed another flaw, namely that it was integrating the negative part of that same function as well as the positive. I had forgotten to implement that Maple should only integrate the positive part of the function. 

So I thought the easiest way for Maple to integrate the function would be to use piecewise() to get rid of the negative part and then integrate the piecewise function. To test this, I temporarily changed the integration interval so that no discontinuity would be present. However, Maple runs for an infinite amount of time when it tries to integrate the piecewise function. Also, I still don't have a clue how to deal with the discontinuities that occur in the real integration interval (when using piecewise). 

Attached at the bottom of this message is the code taken from the inner most for loop. In italic it shows the "problem section" with the piecewise function. The plots were there to help me visualise the functions and what's happening. They'll be removed in the final code.

Can someone help me to try and get the integration of the piecewise function to work? The integration interval in the code is the real integration interval. The discontinuities occur there where peicewise makes the function zero.

Many thanks!

restart: with(plots): with(linalg): with(Optimization): with(ListTools):

Average percentage of sunhours per day in a certain month
Sun_jan:=0.735: Sun_feb:=0.715: Sun_mar:=0.67:
Sun_apr:=0.615: Sun_may:=0.64: Sun_jun:=0.56:
Sun_jul:=0.43: Sun_aug:=0.37: Sun_sep:=0.43:
Sun_oct:=0.715: Sun_nov:=0.81: Sun_dec:=0.74:

Constants
I_sc:=1367.7:
H:= 150:
t_D:=0:
lat:=evalf((9+24/60+27/3600)*Pi/180):
lon:=evalf((51/60+12/3600)*Pi/180):
LC:=evalf((lon*180/Pi-0)/15):
alfa_alt:=evalf(-0.83*Pi/180):

H_zeni:=1.2:
H_azi:=3.6:
N:=116:

E_year:=0:

Initialise all necessary collector and radiation variables

I_o:=evalf(I_sc*(1+0.0333*cos((2*Pi*N)/365))):
delta:=evalf(arcsin(0.39795*cos(0.98563*(N-173)*Pi/180))):
x_input:=(2*Pi*(N-1))/365.242:
EOT:=x->evalf(0.258*cos(x)-7.416*sin(x)-3.648*cos(2*x)-9.228*sin(2*x)):
ts_input:=t->t+EOT(x_input)/60-LC-t_D:
omega:=ts->15*(ts-12)*Pi/180:
omega_set:=evalf((arccos((sin(alfa_alt)-sin(lat)*sin(delta))/(cos(lat)*cos(delta))))):
ts_rise:=evalf(-(omega_set*180/Pi)/15+12):
ts_set:=evalf((omega_set*180/Pi)/15+12):
t_rise:=ts_rise-EOT(x_input)/60+LC+t_D:
t_set:=ts_set-EOT(x_input)/60+LC+t_D:
t_noon:=12-EOT(x_input)/60+LC+t_D:
alfa_elev:=omega->evalf(arcsin(sin(delta)*sin(lat)+cos(delta)*cos(lat)*cos(omega))):
theta:=alfa_elev->Pi/2-alfa_elev:
alfa_azitest:=evalf(arccos((sin(delta)*cos(lat)-cos(delta)*cos(omega(ts_input(t)))*sin(lat))/cos(alfa_elev(omega(ts_input(t)))))):

Calculate direct and normal intensity I with mathematical model on time t of day N

a_0:=0.4237-0.00821*(6-H*0.001)^2:
a_1:=0.5055+0.00595*(6.5-H*0.001)^2:
k_HCDmod:=0.2711+0.01858*(2.5-H*0.001)^2:
I_HCDmod:=theta->I_o*(0.95*a_0+0.98*a_1*exp(-(1.02*k_HCDmod)/cos(theta))):

theta_noon:=evalf(theta(alfa_elev(0))):
I_noon:=evalf(I_HCDmod(theta_noon))*3600: Zet I om van W/m² naar J/(m².uur) voor toekomstige integratie
I_HSmod:=evalf(I_noon*sin(((t-t_rise)*Pi)/(t_set-t_rise)));

Calculate the real radiation energy on the collector for a whole clear day N per (m²)

RF_morning:=evalf(sin(alfa_elev(omega(ts_input(t))))*cos(H_zeni)+cos(alfa_elev(omega(ts_input(t))))*sin(H_zeni)*cos(H_azi-alfa_azitest));
plot(RF_morning,t=t_rise-1..t_noon);

RF_afternoon:=evalf((sin(alfa_elev(omega(ts_input(t))))*cos(H_zeni)+cos(alfa_elev(omega(ts_input(t))))*sin(H_zeni)*cos(H_azi-(2*Pi-alfa_azitest))));
plot(RF_afternoon(t),t=t_noon..t_set+1);

I_real:=piecewise(RF_morning>=0 and I_HSmod>=0 and t<=t_noon,RF_morning*I_HSmod,RF_afternoon>=0 and I_HSmod>=0 and t>t_noon,RF_afternoon*I_HSmod);
plot(I_real,t=t_rise-1..t_set+1);

E_clear:=int(I_real,t=t_rise..t_set,numeric=true);

Apply the average sunhours per day to the radiation energy to go get the real I for an average day (same percentage per hour)

if N<=31 then E_av:=E_clear*Sun_jan elif (N>31 and N <=59) then E_av:=E_clear*Sun_feb
elif (N>59 and N <=90) then E_av:=E_clear*Sun_mar elif (N>90 and N <=120) then E_av:=E_clear*Sun_apr
elif (N>120 and N <=151) then E_av:=E_clear*Sun_may elif (N>151 and N <=181) then E_av:=E_clear*Sun_jun
elif (N>181 and N <=212) then E_av:=E_clear*Sun_jul elif (N>212 and N <=243) then E_av:=E_clear*Sun_aug
elif (N>243 and N <=273) then E_av:=E_clear*Sun_sep elif (N>273 and N <=304) then E_av:=E_clear*Sun_oct
elif (N>304 and N <=334) then E_av:=E_clear*Sun_nov elif (N>334 and N <=365) then E_av:=E_clear*Sun_dec end if:

E_av;