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)?

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

Answer 1

[Walter Roberson]

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.

[mustapha mustapha]

thank you Walter for your answer

I understood from your answer that I should use event functions instead of using if.

I will try that , thanks (:

Answer 2

[Alan Stevens]

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.

[mustapha mustapha]

thank you Alan for your answer

it seems I will use event functions as Walter said.

because the third condition can not be taken out (the temperature has a cyclical behavior)