how to change the system of ODEs (using event functions) when one of the solutions reaches certain value?

Question

mustapha mustapha 2021-1-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/706513-how-to-change-the-system-of-odes-using-event-functions-when-one-of-the-solutions-reaches-certain-v

评论： mustapha mustapha 2021-1-5

Hello

I have a system of 10 coupled differential equations (with initial conditions), and in this system there are parameters depending on the solutions of this system, and other parameters are constant. the thing I want to do is to change the system when one of the solutions (y (2)) reaches a certain value. so the final values of the solutions of the first system will be the initial conditions of the following system .

I used the ode45 solver. I have used event functions, but changing the system doesn't work. may be I didn't use it correctly.

I want the first change occurs when y(2)>=638.9 , and the second change occurs when y(2)<638.9

thank you in advance.

this is the first system (the difference between systems is in dydt(1), dydt(2) and dydt(3)):

%===================param depending on y( )==============
rho_w = (-3.51e-6 * y(2)^2 + 2.01e-3 * y(2) + 0.7125);
hC = 4.18 * rho_w * y(1)+ 1.8364e54;
Pw = y(4) * X_f / 1e-4;
rho = (1/beta) - ((1/beta) * (eta * 2.3102e-23 / 2.3102e-23 + gamma_mf * ((0.1111*y(1)*rho_w) / 2.8867e+07) ... 
            * exp(-gamma_mp / ( 2.8867e+07 / (0.1111*y(1)*rho_w) )) ...
           * epsilon * P_F_fast * P_F_thermal + gamma_w * y(1) / V_total));
%================= dif eqts===========================
dydt(1) = 1.5779e+10*(1.8607e+07-y(1))/(1.2355e+07);
dydt(2) = (1 / hC ) * (Pw - (0.0095 * (y(2) - 423) * 7.1060e+03 ))...
        - 1.5779e+10 * ( rho_w * 4.18 * ( y(2) - 423) / hC)* ((1.8607e+07 - y(1) ) / (1.2355e+07));
dydt(3) = 0;
dydt(4) = ( beta * y(4) / 1e-4 ) * ( rho - 1 ) ...
    +(  lambda_1 * y(5) ...
      + lambda_2 * y(6) ...
      + lambda_3 * y(7) ...
      + lambda_4 * y(8) ...
      + lambda_5 * y(9) ...
      + lambda_6 * y(10) ) + S_spont ; 
dydt(5) = (beta_1 * y(4)/1e-4) - lambda_1 * y(5);
dydt(6) = (beta_2 * y(4)/1e-4) - lambda_2 * y(6);
dydt(7) = (beta_3 * y(4)/1e-4) - lambda_3 * y(7);
dydt(8) = (beta_4 * y(4)/1e-4) - lambda_4 * y(8);
dydt(9) = (beta_5 * y(4)/1e-4) - lambda_5 * y(9);
dydt(10) = (beta_6 * y(4)/1e-4) - lambda_6 * y(10);
%========================= initial conditions =====================
V0_w = 6.2518e+06;
T0 = 423;
d0_gw = 0;
n_0 = 0;
C0_1 = 0;
C0_2 = 0;
C0_3 = 0;
C0_4 = 0;
C0_5 = 0;
C0_6 = 0;
%
y0 = [V0_w;T0;d0_gw;n_0;C0_1;C0_2;C0_3;C0_4;C0_5;C0_6];