Search code examples
matlabsimulink

MATLAB event location to find time required for full conversion of different species


I have a system of ODEs as follows:

dx1/dt = (x1,x2,x3)
dx2/dt = (x1,x2,x3)
dx3/dt = (x1,x2,x3)

The initial conditions are x1=x2=x3=0 @ t=0 and the constraints are dx1/dt = 0, dx2/dt = 0, dx3/dt = 0 for x1 = 1, x2 = 1, x3 = 1 respectively. Once x1, x2, x3 reach value of 1 they should remain at 1 for further increase in t.

I need to find out (1) x1, x2, x3 at different t and (2) estimate the values of t when each of them becomes 1. I had difficulty in getting ressults for (2).

I tried to use the following event function:

function [val,stop,dir] = event(t,X)
X1 = x1; 
X2 = x2; 
X3 = x3;
val = [X1 -1; X2 -1; X3 -1];
stop= [1;1;1];
dir = [1;1;1];
end

It did not work. Then I tried another code to at least find t corresponding to x3 = 1 (because x3 increases slowly compared to x1 and x2).

function [val,stop,dir] = event(t,X)
X3 = x3;
val = X3 -1;
stop= 1;
dir = 1;
end

Could anyone guide me in this regard?

With regards. Sudip


Solution

  • This looks like pretty straightforward to solve with the ode solvers, without the need for event. I would define a function such as:

    function dX = my_ode(t,X)
    
    for k=1:length(X)
        if abs(X(k)-1) <= 1e-6 % or use whatever tolerance you want, e.g. eps
           dX(k) = 0;
        else
           dX(k) = X(k);
        end
    end
    

    I would then call the ode solver as follows:

    tspan = [0 10]; % choose whatever time span you want
    X_init = [0 0 0];
    [T,X] = ode45(@my_ode,tspan,X_init);
    plot(T,X)
    

    Not sure why you have the Simulink tag, but that's also easily implemented in Simulink if need be.