MATLAB: I am trying to solve ODE having more than 1 dependent variable. but I am not able to solve it with dsolve function. please help me to find correct function to solve ODE.This is an mechanical engineering equation so it is bit complicated and lengthy

dsolveMATLABmatlab coderMATLAB Online Serverodeode113ode23ode45

syms Temp_w Temp_g Temp_col
%defining constant
m_w = 15; % mass of water in Kg
m_g = 5; % mass of glass cover in Kg
m_col = 6; % mass of collecter in Kg
c_w = 4179.6; % specific heat of water
c_g = 800; % specific heat of glass
c_col = 2000; % specific heat of collecter
a_w = 0.05; % absorbtivity of water
a_g = 0.08; % absorbtivity of glass
a_col = 0.9; % absorbtivity of collecter
Temp_a = 300; % ambient temerature in Kelvin
G_max = 630; % maximum irradiation of sun
G = 1:10;
for time = 1:10
G(time) = (G_max/2)*(sin((pi*time)/12));
end
G_avg = mean(G); % avg of solar irradiation during a day
% belove are term use in ODE equation
q_ev = 16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))));
q_cw = 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g);
q_rw = 5.67 * (10^(-8)) * 0.96 * ((Temp_w^4) - (Temp_g^4));
q_ca = 6.4 * (Temp_g - Temp_a);
q_ra = 5.67 * (10^(-8)) * 0.93 * ((Temp_g^4) - (Temp_a^4));
q_w = 1000 * (Temp_col - Temp_w);
q_ins = (0.026 * (Temp_col - Temp_a))/0.01;
ag_G = a_g*G_avg;
aw_G = a_w*G_avg;
acol_G = a_col*G_avg;
% ODE equation(1) = q_ev + q_cw + q_rw + ag_G - q_ra - q_ca == m_g*c_g* (dTemp_g/dt)
syms Temp_g(t)
ode_1 = m_g*c_g*diff(Temp_g,t) == 16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) + 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g) + 5.67 * (10^(-8)) * 0.96 * ((Temp_w^4) - (Temp_g^4)) + (a_g*G_avg) - 6.4 * (Temp_g - Temp_a) - 5.67 * (10^(-8)) * 0.93 * ((Temp_g^4) - (Temp_a^4));
Temp_g_sol(t) = dsolve(ode_1)
% ODE equation(2) = -q_ev - q_cw - q_rw + aw_G + q_w == m_w*c_w* (dTemp_w/dt)
syms Temp_w(t)
ode_2 = m_w*c_w*diff(Temp_w,t) == -(16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))) - (0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g)) - (5.67 * (10^(-8)) * 0.96 * ((Temp_w^4) - (Temp_g^4))) + (a_w*G_avg) + (1000 * (Temp_col - Temp_w));
Temp_w_sol(t) = dsolve(ode_2)
% ODE equation(3) = acol_G - q_w - q_ins == m_col*c_col* (dTemp_col/dt)
syms Temp_col(t)
ode_3 = m_col*c_col*diff(Temp_col,t) == (a_col*G_avg) - (1000 * (Temp_col - Temp_w)) - ((0.026 * (Temp_col - Temp_a))/0.01);
Temp_col_sol(t) = dsolve(ode_3)

Best Answer

Here's a possible way using ode45. I've used arbitrary initial values for the three temperatures, so you will need to replace them with the true values.
a_w = 0.05; % absorbtivity of water
a_g = 0.08; % absorbtivity of glass
a_col = 0.9; % absorbtivity of collecter
G_max = 630; % maximum irradiation of sun
G = 1:10;
for time = 1:10
G(time) = (G_max/2)*(sin((pi*time)/12));
end
G_avg = mean(G); % avg of solar irradiation during a day
ag_G = a_g*G_avg;
aw_G = a_w*G_avg;
acol_G = a_col*G_avg;
absorp = [ag_G, aw_G, acol_G];
tspan = [1 10];
Temp_g0 = 310; %%% ? Replace with true values.
Temp_w0 = 320; %%% ?

Temp_col0 = 330; %%% ?
% Combine the three initial temperatures into a single vector.
T0 = [Temp_g0, Temp_w0, Temp_col0];
% Call the function that calculates the rates of change with ode45
% using the following syntax (try doc ode45 for more details).
[t, T] = ode45(@(t,T) fn(t,T,absorp),tspan,T0);
% Extract the individual values of the three temperatures from T.
Temp_g = T(:,1); Temp_w = T(:,2); Temp_col = T(:,3);
% Display the results
plot(t,Temp_g,t,Temp_w,t,Temp_col),grid
xlabel('time'),ylabel('Temps')
legend('Tg','Tw','Tcol')
% The rate of change of temperatures wrt time are contained in the following function.
% The input arguments must be t and T or ~ and T (the latter if t is not used explicitly
% within the function). Also, any extra data that is used within the function
% ca be passed in as an extra argument - like absorp here.
% The three temperatures are passed in via a single vector T.
% Their individual rates of change are returned as a column vector in dTdt.
function dTdt = fn(~,T,absorp)
% Extract the three individual temperatures for use in the rates
% of change equations below.
Temp_g = T(1); Temp_w = T(2); Temp_col = T(3);
% Extract the individual pieces of data passed to the function
% in absorp.
ag_G = absorp(1);
aw_G = absorp(2);
acol_G = absorp(3);
%defining constant
m_w = 15; % mass of water in Kg
m_g = 5; % mass of glass cover in Kg
m_col = 6; % mass of collecter in Kg
c_w = 4179.6; % specific heat of water
c_g = 800; % specific heat of glass
c_col = 2000; % specific heat of collecter
Temp_a = 300; % ambient temerature in Kelvin
% ODE equation(1) = q_ev + q_cw + q_rw + ag_G - q_ra - q_ca == m_g*c_g* (dTemp_g/dt)
dTemp_gdt = (16.273 * 0.884 * (((Temp_w - Temp_g + ...
(((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - ...
(3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - ...
(3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - ...
(3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) *...
(11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) +...
0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - ...
(11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - ...
(3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - ...
(11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g) + 5.67 * (10^(-8)) * ...
0.96 * ((Temp_w^4) - (Temp_g^4)) + ag_G - 6.4 * (Temp_g - Temp_a) - ...
5.67 * (10^(-8)) * 0.93 * ((Temp_g^4) - (Temp_a^4)))/(m_g*c_g);
% ODE equation(2) = -q_ev - q_cw - q_rw + aw_G + q_w == m_w*c_w* (dTemp_w/dt)
dTemp_wdt = (-(16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - ...
(3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*...
(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) *...
(11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) *...
(11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))) - ...
(0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - ...
(11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - ...
(3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - ...
(11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g)) - (5.67 * (10^(-8)) * ...
0.96 * ((Temp_w^4) - (Temp_g^4))) + aw_G + (1000 * (Temp_col - Temp_w)))/(m_w*c_w);
% ODE equation(3) = acol_G - q_w - q_ins == m_col*c_col* (dTemp_col/dt)
dTemp_coldt = (acol_G - (1000 * (Temp_col - Temp_w)) - ...
((0.026 * (Temp_col - Temp_a))/0.01))/(m_col*c_col);
% Assemble a column vector of the values of the rates of change.
dTdt = [dTemp_gdt; dTemp_wdt; dTemp_coldt];
end