Confusion on event function

Hi all, I want to incorporate an event into my ode45 code. The event is that every time the value y reaches 1 (initial condition 0) it should stop integrating and must restart from zero. I will share the code I have written below. Looks fine, but I am not able to understand how it works. Actually where did I set the event in this code? (I was just imitating a code given in the documentation). If it works for N=1, I need to generalise it for large number of systems. Will that be possible?

N=1;                                                                     
m=0.3;                                          

   %Initial values                                                                      
  ystart=zeros(N,1);

options = odeset('Events',@SpikeEvents,'RelTol',1e-5,'AbsTol',1e-4);

%ODE solver
[t,y,te,ye,ie]=ode45(@ELIF,tspan,ystart,options);

%Plotting
plot(t,y)
xlabel('time')
ylabel('voltage')
title('Leaky IF with spike reset')

%Defining event

%Defining event
function [value,isterminal,direction] = SpikeEvents(t,y)
value = y(1);
isterminal = 1;   % Stop the integration
direction = 1;   
end

%Defining the system of equations

function dydt=ELIF(t,y,m)
N=1;
m=0.3;
dydt=zeros(N,1);
for i=1:N
    dydt(i)=-y(i)+m; %The model equation       

end
end

Answer 1

[Walter Roberson]

Event functions cannot reset the integration by themselves. Instead what you need to do is configure the event function to signal termination, and then you need to loop, something like

K = 1;
start_time = tspan(1);
end_time = tspan(2);
time_tolerance = 1e-3;
y0 = ystart;

while true
  [t{K}, y{K}, te{K}, ye{K}, ie{K}]=ode45(@ELIF, [start_time end_time], ystart, options);
  ended_at = t{K}(end);
  if end_time - ended_at <= time_tolerance
    break;
  end
  start_time = ended_at;
  y0 = y{K}(end,:);
  y0(1) = 0;   %restart integration from 0
  K = K + 1;
end

Then the overall path would be found be concatenating together all of the t{:} and y{:}, except omitting the overlap time point.

[Vipin Padinjarath]

Thank you Walter. I will get back once I work it out.

[Vipin Padinjarath]

I understand the logic. Now my confusion is, in the above code, where the the terminatikn condition specified? That is I need to stop integration everytime the numerical value of the dependent variable reaches 1.

[Torsten]

%Defining event
function [value,isterminal,direction] = SpikeEvents(t,y)
  value = y(1)-1.0;
  isterminal = 1;   % Stop the integration
  direction = 1;   
end

[Vipin Padinjarath]

Thanks again?

[Steven Lord]

The standard simple example I point people to when trying to explain ODE solver events is the ballode example. It's similar to the one Walter posted, though it grows the arrays of output in its for loop rather than using cell arrays.

[Vipin Padinjarath]

I have gone through the ballode code. Yet I could not fully understand how it works. Thanks for the input. These are great helps.

Answer 2

[Vipin Padinjarath]

In Walter's code, we need to guess the end time? I mean do we need to guess the time when the first event would occur?

[Walter Roberson]

No, just use the last element of your current tspan variable, the point where you would want to end anyhow.

[Vipin Padinjarath]

Ok. ?

[Vipin Padinjarath]

Thank you once again Walter. Finally it worked the way I wanted it to work. The looping works perfectly fine.

[Vipin Padinjarath]

Hi, I wanted to extend the above code for a network of N such systems. I want all all the component of the system to reset to zero when y reaches 1 and the starts integrating again. The code runs; but with following problems;

1. I am running it for N such systems. It looks as if the integration stops when the first system (N=1, y=y(1)) meets the event and restarts. It doesn't check whether the other systems have crossed the event. Because in the numerical values it generates, though first system never crosses 1, others are crossing 1.
2. Then I tried to couple them through simple means. The coupling doesn't seem to have any impact on the code. The numerical values are expected to change when coupling is introduced. Once the resetting is done and starts the integration again from zero, all the systems evolves identically; doesn't matter whether there is xoupling or not.

I have beein trying to fix this for las 40-45 days. The code is attached here. Your help have been great and I owe you alot for helping me learning matlab better.

                             
N=10; %Number of neurons
tspan=(0:10);   %time steps                                                                                                       
tstart=tspan(1); %first element of tspan vector as the starting time of integration
tend=tspan(end);
y0=rand(N,1); %random initial values of neurons

options = odeset('Events',@SpikeEvents,'RelTol',1e-5,'AbsTol',1e-4);

%Inidtializing cells to store values
k=1;
t{k}=[];
y{k}=[];
te{k}=[];
ye{k}=[];
ie{k}=[];

%solving the equation
while k<length(tspan)
    [t{k}, y{k}, te{k}, ye{k}, ie{k}]=ode45(@CEnIF,[tstart tend], y0, options);
    endtime=t{k}(end); %end time is the last element of the kth cell of time vector. k th cells contains time steps from tstart to te(time step where event occures)
    tstart=endtime+0.002;%resetting happens after a small time interval
    y0=y{k}(end);
    y0=zeros(N,1); %restart the integration with new initial value
    k=k+1;
end

%Defining event

function [value,isterminal,direction] = SpikeEvents(~,y)
N=10;
for i=1:N
value = y(1)-1;
isterminal = 1;   % Stop the integration
direction = [];   
end
end

%Defining the system of equations

function dydt=CEnIF(~,y,m,N,S)
N=10;
m=0.99;
S=0.5; %coupling streangth
dydt=zeros(N,1);
for i=1:N
    if i<N
        dydt(i)=exp(y(i))+m+(S/6)*sum((y(i)-y(i+1)));
    else
       break
    end
end    
end
   

    

[Walter Roberson]

why loop to N and have the if? Why not loop to N-1 ?

You would not even need to loop just vectorize.

[Walter Roberson]

value = y - 1;

would appear to provide testing for all of the systems .

[Vipin Padinjarath]

I replaced the loop within the event function to y(1)-1. Still one of the yis crossing the event.

[Vipin Padinjarath]

looping to N-1 is also done. But as I said the coupling term has no impact on the system itself.

[Walter Roberson]

no loop just

value = y - 1;

[Vipin Padinjarath]

y-1 or y(1)-1?

[Walter Roberson]

y-1

if you want to test all of the y values at the same time which is what I have interpreted your question about simultaneous systems to be .

[Vipin Padinjarath]

Okay. Thank you very much.