1 个投票

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

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2018-10-26

1 个投票

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.

6 个评论
显示 4更早的评论 隐藏 4更早的评论

Vipin Padinjarath
Vipin Padinjarath 2018-10-26
Thank you Walter. I will get back once I work it out.
Vipin Padinjarath
Vipin Padinjarath 2018-10-26
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
Torsten 2018-10-26
%Defining event
function [value,isterminal,direction] = SpikeEvents(t,y)
value = y(1)-1.0;
isterminal = 1; % Stop the integration
direction = 1;
end
Steven Lord
Steven Lord 2018-10-26
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
Vipin Padinjarath 2018-10-26
I have gone through the ballode code. Yet I could not fully understand how it works. Thanks for the input. These are great helps.

请先登录，再进行评论。

更多回答（1 个）

Vipin Padinjarath
Vipin Padinjarath 2018-10-30

1 个投票

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?

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Walter Roberson
Walter Roberson 2018-10-30
No, just use the last element of your current tspan variable, the point where you would want to end anyhow.
Vipin Padinjarath
Vipin Padinjarath 2018-11-14
Thank you once again Walter. Finally it worked the way I wanted it to work. The looping works perfectly fine.
Vipin Padinjarath
Vipin Padinjarath 2018-12-3
编辑：Vipin Padinjarath 2018-12-3
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
Walter Roberson 2018-12-3
编辑：Walter Roberson 2018-12-3
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
Walter Roberson 2018-12-3
value = y - 1;
would appear to provide testing for all of the systems .
Vipin Padinjarath
Vipin Padinjarath 2018-12-3
I replaced the loop within the event function to y(1)-1. Still one of the yis crossing the event.
Vipin Padinjarath
Vipin Padinjarath 2018-12-3
looping to N-1 is also done. But as I said the coupling term has no impact on the system itself.
Walter Roberson
Walter Roberson 2018-12-3
no loop just
value = y - 1;
Walter Roberson
Walter Roberson 2018-12-8
编辑：Walter Roberson 2018-12-12
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
Vipin Padinjarath 2018-12-12
Okay. Thank you very much.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by