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
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
RUTVIK RAVAL
el 23 de Feb. de 2021
Comentada: Alan Stevens
el 24 de Feb. de 2021
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)
3 comentarios
Respuesta aceptada
Alan Stevens
el 23 de Feb. de 2021
Editada: Alan Stevens
el 24 de Feb. de 2021
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
2 comentarios
Alan Stevens
el 24 de Feb. de 2021
Ok. I've added a few comments. You should also study the documentation on ode45.
Más respuestas (0)
Ver también
Categorías
Más información sobre Symbolic Math Toolbox en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!