HOW TO SOLVE A SYSTEM OF FIRST ORDER ODE WITH ODE45 SOLVER (knowing that it contains if statment and other parameter depending on the solutions of the system)?
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
mustapha mustapha
el 23 de Dic. de 2020
Comentada: mustapha mustapha
el 23 de Dic. de 2020
HELLO
I would I would be grateful if you could help me.
i have a system of ODE that i want to solve with ode45;
i used if statment but it didn't seem to work properly (it gives the solution for just one condition)
the function is : function F = ODEfun(t,y)
the system is given below
i defined the ode45 solver:
[t, y] = ode45(@ODEfun,tspan,y0);
F = [dydt(1);dydt(2);dydt(3)];
and this is the system
%===========================================================================================
rho_w = a_1 * y(2)^2 + b_1 * y(2) + c_1;
%=========================================================================================
hC = a_2*y(1) + b_2 + c_2 + d_2;
% ================================ Differential equations =========================================
% ================================ water content =========================================
if y(2)<T_b
dydt(1) = Q*(phi*V_total-y(1))/(phi_eff*V_total);
elseif y(2)>= T_b
dydt(1) = -(h*A_s/(rho_w * L_v))*(y(2) - T_b);
elseif y(2)<T_b
dydt(1) = 0;
end
%========================================== temperature===================
if y(2) < T_b
dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl / d_inf ))...
- Q * ( rho_w * c_w * ( y(2) - T_inf)/hC)* ((phi * V_total - y(1) ) /...
phi_eff*V_total);
elseif y(2)>= T_b
dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl / d_inf ))...
- ( rho_w * L_v * h * A_cyl * (y(2) - T_b) / (hC * rho_w * L_v));
elseif y(2)<T_b
dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl / d_inf ));
end
%=============================================Distance ======================
if y(2)<T_b
dydt(3) = 0;
elseif y(2)>= T_b
dydt(3) = (h * A_s / (rho_w * L_v * A ) ) * ( y(2) - T_b );
elseif y(2)<T_b
dydt(3) = -Q / A;
end
0 comentarios
Respuesta aceptada
Walter Roberson
el 23 de Dic. de 2020
Alan Stevens pointed out the conflict in your conditions. However, the only time that y(2) < T_b and y(2)>= T_b can both be false is if y(2) is nan, so your third condition would not have been reached anyhow.
However:
The mathematics of Runge Kutta assumes that the derivatives you give have continuous first and second derivatives. When you use if inside your functions, it is seldom the case that the first and second derivatives of each returned value is continuous. You would have had to have arranged your
if y(2)<T_b
dydt(1) = Q*(phi*V_total-y(1))/(phi_eff*V_total);
elseif y(2)>= T_b
dydt(1) = -(h*A_s/(rho_w * L_v))*(y(2) - T_b);
elseif y(2)<T_b
dydt(1) = 0;
end
to have matched the first and second derivatives of Q*(phi*V_total-y(1))/(phi_eff*V_total) and -(h*A_s/(rho_w * L_v))*(y(2) - T_b) at T == T_b, and match 0 at whatever the third condition is.
It isn't impossible... it is just messy. It can happen in some cases involving cubic spline interpolation though.
If you cannot match two derivatives your expression branches at the condition breakpoints, then ode45() will not be able to calculate the ODE properly. It will either notice the problem and complain about a singularlity, or else it will not notice the problem and will silently return the wrong value, which is a much worse circumstance (because it misleads you into thinking your answer is not nonsense when it is nonsense.)
If you cannot match the derivatives then you need to use event functions that signal termination when the conditions are crossed, and then you have to restart integration from where you left off.. this time on the other side of the boundary.
Más respuestas (1)
Alan Stevens
el 23 de Dic. de 2020
Look at the two asterisked lines in your temperature section below
if y(2) < T_b %********************************************************************
dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl / d_inf ))...
- Q * ( rho_w * c_w * ( y(2) - T_inf)/hC)* ((phi * V_total - y(1)) /...
phi_eff*V_total);
elseif y(2)>= T_b
dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl / d_inf ))...
- ( rho_w * L_v * h * A_cyl * (y(2) - T_b) / (hC * rho_w * L_v));
elseif y(2)<T_b %********************************************************************
dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl / d_inf ));
end
The conditions are identical, so the last elseif will never be implemented!
The same is true of your water content and distance sections.
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!