How to insert an if conditional in ode45?

Hi, I'm trying to simulate a leaky integrate-and-fire model. I'm using Ode45, it all works well but I want to insert an If conditional so that when the y value is less than a specific value (-50), y must reset to the value -65. I have no idea about how to implement this. Any suggestions? Thanks a lot.

x0=-70;
t0=0;
dt=0.1;
tmax=2;
tvec=t0:dt:tmax;
[t,x]=ode45(@test_ode, tvec, x0)
plot (t,x)

function [dydt]=test_ode(t,y)
         dydt=-1/10*y+13     
end 

Answer 1

[Star Strider]

Try this:

x0=-70;
t0=0;
dt=0.1;
tmax=20;
% tvec=t0:dt:tmax;
tvec = [t0 tmax];
opts = odeset('Events',@depolEvents);

for k = 1:10
    [t{k},x{k}]=ode45(@test_ode, tvec, x0, opts);
    x0 = -65;
    tvec = [t{k}(end) t{k}(end)+tmax]
%     tvec = t{k}(end) : dt : t{k}(end)+tmax;
end

tv = cat(1,t{:});
xv = cat(1,x{:});

figure
plot (tv,xv)
grid

function [dydt]=test_ode(t,y)
         dydt=-1/10*y+13;     
end  

function [position,isterminal,direction] = depolEvents(t,y)
    position = y(1) + 50; % The value that we want to be zero
    isterminal = 1;  % Halt integration 
    direction = 1;   % The zero can be approached from either direction
end

This uses the techniques in ODE Event Location to stop the integration that is restartted in the next increment of the loop.

If it does not do what you want as written (since I may not have understood exactly what you want to do), experiment with it to get the result you want.

[Simone Maucci]

Thank you very much for your answer. I tried it but unfortunately I got the following error: "Unable to perform assignment because the left and right sides have a different number of elements", concerning line 10:

    [t{k},x{k}]=ode45(@test_ode, tvec, x0, opts);

Since I'm not familiar with this function, can you help me to solve it? thanks a lot

[Star Strider]

My pleasure!

My posted code ran without error in R2020b. I tested it before I posted it.

I have no idea what could be causing that problem, since the ‘left sides’ are cell arrays, and so can be anything, or vectors or matrices of any size. The components of ‘t’ and ‘x’ are (or should be) column vectors in this code, and were when I tested it.

When I ran it using the posted code, it produced this plot:

EDIT — (1 Feb 2021 at 19:40)

I still do not understand the problems you are having with the cell arrays.

This is a work-around that avoids them:

tv = [];
xv = [];
for k = 1:10
    [t,x]=ode45(@test_ode, tvec, x0, opts);
    x0 = -65;
    tvec = [t(end) t(end)+tmax]
    tv = [tv; t];
    xv = [xv; x];
end

figure
plot (tv,xv)
grid

Use that instead, and delete this complete block of code:

for k = 1:10
    [t{k},x{k}]=ode45(@test_ode, tvec, x0, opts);
    x0 = -65;
    tvec = [t{k}(end) t{k}(end)+tmax]
end
tv = cat(1,t{:});
xv = cat(1,x{:});

Since a partial replacement will throw errors, delete the block using the cell arrays and replace it with the first block in this edited Comment.

The full revised code is now:

x0=-70;
t0=0;
dt=0.1;
tmax=20;
tvec = [t0 tmax];
opts = odeset('Events',@depolEvents);

tv = [];
xv = [];
for k = 1:10
    [t,x]=ode45(@test_ode, tvec, x0, opts);
    x0 = -65;
    tvec = [t(end) t(end)+tmax]
    tv = [tv; t];
    xv = [xv; x];
end

figure
plot (tv,xv)
grid

function [dydt]=test_ode(t,y)
         dydt=-1/10*y+13;     
end  

function [position,isterminal,direction] = depolEvents(t,y)
    position = y(1) + 50; % The value that we want to be zero
    isterminal = 1;  % Halt integration 
    direction = 1;   % The zero can be approached from either direction
end

That works correctly, and produces the same figure.

.

[Simone Maucci]

Now it works perfectly, no sign of errors anymore. I'll try to practice on it a little bit. Thank you very much for your kindness, I really appreciated it.

[Star Strider]

As always, my pleasure!

Answer 2

[Steven Lord]

Don't use a conditional. Use an events function to detect when the condition occurs, stop the ODE solver, use the result from the time of the event to create initial conditions for a new call to the solver, and call the ODE solver from the time of the event.

You should be able to use the ballode example as a model for your function, though instead of bouncing a ball off the floor at ground level you're "bouncing it off the floor" at -50.

[Simon Mwakitabu]

Now imagine you have got a condition which depends on the variable itself. For example you have got a check valve equation, where it opens when P1(t)*Area-P2(t)*Area>k*x_o While solving lets say dP1dt=Q_check.... dP2dt=-Q_check....

[Steven Lord]

The events function in the ballode example detects when the height of the ball above the ground (which is the variable for which we're solving) crosses 0 (representing the ball hitting the ground.)

If you're having difficulty writing up an events function for your problem you could ask a new question (linking back here for reference) showing the code you've written to try to solve the problem as well as the mathematical form of the differential equations you're trying to solve and the posters on Answers may be able to give you some advice.