Is it possible to have an ode solver solve in chronological order?

I'm trying to solve a piecewise differential equation of Temperature as a function of time using ode45. Some of the conditions of the piecewise are dependent on timestamps of certain events in the differential function. For now, I am using global variables to store those times. The plot did not behave the way I expected, so I displayed time in the ode function, and saw that time was not monotonically increasing. This caused the timestamps to be changed after calculations dependent on them were performed. I looked into event location functions, but if ode45 does not solve in order, I don't think that'll work either. Is it possible to use an ode solver for this problem?

Here's an example of the script:

f = @(t,T)funcOfT(t,T)
[t,T] = ode45(f, [0 6000], 0);
plot (t,T)

function dTdt = funcOfT(t,T)
    global t1 t2;

    T_max = 5;
    x = 1;
    a = 1;
    b = 1;

    %   Displaying time for troubleshooting
    disp(t);

    if(t<500)
        %   First case ends at timestamp t1
        if (T < T_max)
            dT = T/2;
            t1 = t;
        %   Third case begins at t2
        elseif (t-t1 > x)
            dT = T/3;
        %   Second case starts at t1, ends at t2
        else
            dT = 0;
            t2 = t;
        end
    else
        %   Fourth case ends at t3
        if (T > T_max)
            dT = T/2;
            t3 = t;
        %   Sixth case begins with this timestamp-based condition
        elseif ((t-t3)*a + (t2-t1)*b < 0)
            dT = T/4;
        %   Fifth case begins at t3
        else    
            dT = 0;
        end
    end

end

[Steven Lord]

Can you show us the mathematical form of the differential equations you're trying to solve? Use the third icon (Σ) in the Insert section of the Toolstrip to enter them using TeX or use the first icon in that section to include a picture of the equations. Seeing the mathematical equations (rather than your specific implementation) may allow us to suggest a different approach that is not apparent from your implementation.

[Star Strider]

I would create the timestamps as a double vector (perhaps using datenum) and pass that as the ‘tspan’ vector to the differential equation integrator, such as ode45. I would avoid using global variables at all costs, and instead pass those data as extra parameters.

[Isabella]

For the first, third, fourth, and sixth case, the equation is some form of:

For the second and fifth case, the temperature should be a known constant and

. These two cases are only true if and , respectively.

[Sam Chak]

Hi @Isabella

Can you display the multi-level conditional cases like the following?

It helps people to interpret your problem formulation correctly. I noticed there is another level for cases when the Temperature and .

[Isabella]

For some context, I'm trying to model the temperature of a material as it changes phases and different amounts of power are applied over time. For , total power is positive. During this period, temperature either increases or remains at melting temperature while power input goes towards overcoming latent heat. The rate at which temperature changes depends on whether the temperature should be in the solid or liquid region, which is part of the reason why I want to record and as when the melting region starts and ends. During this latent region, there is no change in temperature, but the product of change in time and the constant power must be tracked. Once this product exceeds the known constant , temperature increases at a different rate.

Then for , the temperature will start either in the liquid region or at melting temperature, depending on the product of difference between and and the constant power applied at that time. Total power is negative for this period. If temperature is at the melting temperature, is set to . Otherwise, temperature will decrease until the it reaches melting temperature. The time at which this happens is tracked as . Now we use to track the amount of heat dissipated, as the difference from current time to mutliplied by the constant power applied must be equal to the heat stored during the period from to in order for the temperature to decrease. Once the heat dissipated exceeds the heat stored, temperature resumes decreasing.

Here's my attempt at writing out the conditionals. (I used generic constants and simplified expressions to make it easier to understand).

Here's where the times come from:

If , .

If , .

If && , .

If , .

Also looking at the ode45 documentation, I do see that it solves using the previous timestep, so I'm now looking into why the displayed times aren't monotonically increasing.

[Sam Chak]

Hi @Isabella

Can you have some realistic values for the parameters? I asked because your example where and initial value which is also the equilibrium point, does not help the simulation.

Request for and the initial temperature , as well as the melting point .

Need to perform analysis on the thermal system.

Since the Temperature (T) is a physical quantity and T > 0 is always true, if and is bounded, then the thermal system is stable.

[Isabella]

Sure, here are some more realistic conditions then:

Initially,

Answer 1

[Bruno Luong]

@Isabella "Also looking at the ode45 documentation, I do see that it solves using the previous timestep, so I'm now looking into why the displayed times aren't monotonically increasing."

If you read Runge Kutta scheme used by ode45 you'll see it evaluate intermediate times in between the so called "time step" then go back the previous time step evalute with the solution of the next time step. So in between two time steps there a few call of the model that are NOT monotonic.

Also ode45 possibly reduces the time step when it estimates the system become stiffer. I don't know the exact detail but this could also means that there is some non monototic calls of your model that jump back the the previus time step.

In your case I would use the event (click on the link) to estimate t1, t2, t3, t4, and NOT detect them inside your model with the assumtion that the model is called with increasing t. That assumption cannot be true.