Why am I not getting desired output?

4 views (last 30 days)
Rakesh
Rakesh on 20 Jun 2022
Commented: Torsten on 20 Jun 2022
I am solving differential equation that contain piecewise function. I coded this in Octave and I am getting desired results but the same in MATLAB is not giving the desired results. Code is compiling without any error but output is different. Code is given below:
function neural_circuit
% using rk4 method
clc;
clear all;
close all;
h = 0.05;
t = 1:h:500;
y0 = [2.5, 3.6, 9.5, 8.8, 2, 3, 7, 7.5];
% t(1) =0;
y = zeros(length(t),8);
y(1,:) = y0;
for i = 1:length(t)
% %updates time
t(i+1) = t(i)+h;
k1 = model(t(i), y(i, :));
k2 = model(t(i)+h/2, y(i, :)+k1*h/2);
k3 = model(t(i)+h/2, y(i, :)+k2*h/2);
k4 = model(t(i)+h, y(i, :)+k3*h);
y(i+1, :) = y(i, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end
plot(t,y(:,1)) %% plot of s_l1
end
function dy= model(t,y)
% dy = zeros(1,8);
%parameters value
tau_s = 10; taur_a = 83; beta_l = 0.0175; Jr_a =148.20; I_ext = 20; I_0=0.29; J_l1l2 = 45.5; J_l1b1= 0.001; J_l2b2 = 0.001; J_b1b2= 44.5;
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;
function y = pieceWise(t)
y = ...
(t ) .* (t >= 0) + ...
(0 ) .* (t <0);
end
% function f = pieceWise(s)
% if s>=0
% f = s;
% else
% f=0;
% end
% end
dy(1) = -(y(1)/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0);
dy(2) = -(y(2)/tau_s)+ beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0);
dy(3) = -(y(3)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0);
dy(4) = -(y(4)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0);
dy(5) = -y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0)./taur_a;
dy(6) = -y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0)./taur_a;
dy(7) = -y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0)./taur_a;
dy(8) = -y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0)./taur_a;
end
  1 Comment
Jan
Jan on 20 Jun 2022
clear all; on top of a function is a complete waste of time and has no meaningful purpose. It deletes all loaded functions from the RAM such that realoading them from the slow disk needs time without any benefit. This is called "cargo cult programming".
"J_l1l1" is really hard to read.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 20 Jun 2022
Edited: Torsten on 20 Jun 2022
h = 0.05;
t = 1:h:500;
y0 = [2.5, 3.6, 9.5, 8.8, 2, 3, 7, 7.5];
% t(1) =0;
y = zeros(length(t),8);
y(1,:) = y0;
for i = 1:length(t)
t(i+1) = t(i)+h;
k1 = model(t(i), y(i, :));
k2 = model(t(i)+h/2, y(i, :)+k1*h/2);
k3 = model(t(i)+h/2, y(i, :)+k2*h/2);
k4 = model(t(i)+h, y(i, :)+k3*h);
y(i+1, :) = y(i, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end
plot(t,y(:,1)) %% plot of s_l1
function dy = model(t,y)
tau_s = 10; taur_a = 83; beta_l = 0.0175; Jr_a =148.20; I_ext = 20; I_0=0.29; J_l1l2 = 45.5; J_l1b1= 0.001; J_l2b2 = 0.001; J_b1b2= 44.5;
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;
dy(1) = -(y(1)/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0);
dy(2) = -(y(2)/tau_s)+ beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0);
dy(3) = -(y(3)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0);
dy(4) = -(y(4)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0);
dy(5) = (-y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0))./taur_a;
dy(6) = (-y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0))./taur_a;
dy(7) = (-y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0))./taur_a;
dy(8) = (-y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0))./taur_a;
end
function y = pieceWise(t)
y = (t ) .* (t >= 0) + ...
(0 ) .* (t <0);
end
  2 Comments
Torsten
Torsten on 20 Jun 2022
dy(5) = -y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0)./taur_a;
dy(6) = -y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0)./taur_a;
dy(7) = -y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0)./taur_a;
dy(8) = -y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0)./taur_a;
instead of
dy(5) = (-y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0))./taur_a;
dy(6) = (-y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0))./taur_a;
dy(7) = (-y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0))./taur_a;
dy(8) = (-y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0))./taur_a;

Sign in to comment.

More Answers (1)

Jan
Jan on 20 Jun 2022
Edited: Jan on 20 Jun 2022
neural_circuit
function neural_circuit
h = 0.05;
t = 1:h:500;
y0 = [2.5, 3.6, 9.5, 8.8, 2, 3, 7, 7.5];
% t(1) =0;
y = zeros(length(t),8);
y(1,:) = y0;
for i = 1:length(t)
t(i+1) = t(i)+h;
k1 = model(t(i), y(i, :));
k2 = model(t(i)+h/2, y(i, :)+k1*h/2);
k3 = model(t(i)+h/2, y(i, :)+k2*h/2);
k4 = model(t(i)+h, y(i, :)+k3*h);
y(i+1, :) = y(i, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end
plot(t,y(:,1)) %% plot of s_l1
end
function dy = model(t,y)
tau_s = 10; taur_a = 83; beta_l = 0.0175; Jr_a =148.20; I_ext = 20; I_0=0.29; J_l1l2 = 45.5; J_l1b1= 0.001; J_l2b2 = 0.001; J_b1b2= 44.5;
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;
function y = pieceWise(t)
y = (t ) .* (t >= 0) + ...
(0 ) .* (t <0);
end
dy(1) = -(y(1)/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0);
dy(2) = -(y(2)/tau_s)+ beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0);
dy(3) = -(y(3)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0);
dy(4) = -(y(4)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0);
dy(5) = -y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0)./taur_a;
dy(6) = -y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0)./taur_a;
dy(7) = -y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0)./taur_a;
dy(8) = -y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0)./taur_a;
end
Output with Octave:
I do not see the problem.
  3 Comments
Rakesh
Rakesh on 20 Jun 2022
You are right. I have coded in octave the same problem using function handle. Please see the octave code below---> this code is giving correct output but I have to code the same by spiliting the whole program into global and local functions that's why i have choosed MATLAB. But in matlab i am not getting the desired results...could you help to split the following program into main and local functions in octave?
% synapses to each sub-population.
graphics_toolkit('gnuplot')
clc;
clear all;
close all;
%defining piecewise functions
function x = pieceWise(s)
x = zeros (size (s));
ind = s > 0;
x(ind) = s(ind);
endfunction
% Defines variables
global tau_s taur_a beta_l I_ext I_0 J_l1l2 J_l1b1 J_l1b2 J_l2b1 J_l2b2 J_b1b2 Jr_a M_l1 M_l2 M_b1 M_b2 J_l1l1 J_l2l2 J_b1b1 J_b2b2;
tau_s = 10; taur_a = 83; beta_l = 0.0175; Jr_a =148.20; I_0=0.29; J_l1l2 = 43.5; J_l1b1= 0; J_l2b2 = 0; J_b1b2= 44.5;
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;
I_ext = 20;
% step sizes
##I_ext = 0:0.001: 20;
t0 = 0;
tfinal = 100;
h = 0.05;
n = ceil((tfinal-t0)/h)+1;
%initial conditions
s_l1(1) = 2.5;
s_l2(1) = 3.6;
s_b1(1) = 9.5;
s_b2(1) = 8.8;
a_l1(1) = 2;
a_l2(1) = 3;
a_b1(1) = 7;
a_b2(1) = 7.5;
t(1) =0;
% functions handles
S_l1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_l1/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*s_l2-J_l1b1*s_b1-J_l1l1*s_l1-a_l1-I_0);
S_l2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_l2/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*s_l1-J_l2b2*s_b2-J_l2l2*s_l2-a_l2-I_0);
S_b1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_b1/tau_s)+ beta_l*pieceWise(I_ext-J_b1b2*s_b2-J_l1b1*s_l1-J_b1b1*s_b1-a_b1-I_0);
S_b2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_b2/tau_s)+ beta_l*pieceWise(I_ext-J_b1b2*s_b1-J_l2b2*s_l2-J_b2b2*s_b2-a_b2-I_0);
A_l1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_l1+Jr_a*beta_l*pieceWise(I_ext-J_l1l2*s_l2-J_l1b1*s_b1-J_l1l1*s_l1-a_l1-I_0))./taur_a;
A_l2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_l2+Jr_a*beta_l*pieceWise(I_ext-J_l1l2*s_l1-J_l2b2*s_b2-J_l2l2*s_l2-a_l2-I_0))./taur_a;
A_b1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_b1+Jr_a*beta_l*pieceWise(I_ext-J_b1b2*s_b2-J_l1b1*s_l1-J_b1b1*s_b1-a_b1-I_0))./taur_a;
A_b2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_b2+Jr_a*beta_l*pieceWise(I_ext-J_b1b2*s_b1-J_l2b2*s_l2-J_b2b2*s_b2-a_b2-I_0))./taur_a;
for i = 1:n
%updates time
t(i+1) = t(i)+h;
%updates S_l1, S_l2, S_b1, S_b2, A_l1, A_l2, A_b1 and A_b2
k1S_l1 = S_l1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1S_l2 = S_l2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1S_b1 = S_b1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1S_b2 = S_b2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1A_l1 = A_l1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1A_l2 = A_l2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1A_b1 = A_b1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
k1A_b2 = A_b2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
% k2th
k2S_l1 = S_l1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2S_l2 = S_l2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2S_b1 = S_b1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2S_b2 = S_b2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2A_l1 = A_l1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2A_l2 = A_l2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2A_b1 = A_b1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
k2A_b2 = A_b2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
% k3
k3S_l1 = S_l1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3S_l2 = S_l2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3S_b1 = S_b1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3S_b2 = S_b2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3A_l1 = A_l1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3A_l2 = A_l2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3A_b1 = A_b1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
k3A_b2 = A_b2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
%k4
k4S_l1 = S_l1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4S_l2 = S_l2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4S_b1 = S_b1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4S_b2 = S_b2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4A_l1 = A_l1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4A_l2 = A_l2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4A_b1 = A_b1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
k4A_b2 = A_b2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
%
s_l1(i+1) = s_l1(i)+h/6*(k1S_l1 + 2*k2S_l1 + 2*k3S_l1 +k4S_l1);
s_l2(i+1) = s_l2(i)+h/6*(k1S_l2 + 2*k2S_l2 + 2*k3S_l2 +k4S_l2);
s_b1(i+1) = s_b1(i)+h/6*(k1S_b1 + 2*k2S_b1 + 2*k3S_b1 +k4S_b1);
s_b2(i+1) = s_b2(i)+h/6*(k1S_b2 + 2*k2S_b2 + 2*k3S_b2 +k4S_b2);
a_l1(i+1) = a_l1(i)+h/6*(k1A_l1 + 2*k2A_l1 + 2*k3A_l1 +k4A_l1);
a_l2(i+1) = a_l2(i)+h/6*(k1A_l2 + 2*k2A_l2 + 2*k3A_l2 +k4A_l2);
a_b1(i+1) = a_b1(i)+h/6*(k1A_b1 + 2*k2A_b1 + 2*k3A_b1 +k4A_b1);
a_b2(i+1) = a_b2(i)+h/6*(k1A_b2 + 2*k2A_b2 + 2*k3A_b2 +k4A_b2);
end
%plots
##%average spike rates
M_l1 = beta_l*pieceWise(I_ext-J_l1l2*s_l2-J_l1b1*s_b1-a_l1-I_0);
M_l2 = beta_l*pieceWise(I_ext-J_l1l2*s_l1-J_l2b2*s_b2-a_l2-I_0);
M_b1 = beta_l*pieceWise(I_ext-J_l1b1*s_l1-J_b1b2*s_b2-a_b1-I_0);
M_b2 = beta_l*pieceWise(I_ext-J_l2b2*s_l2-J_b1b2*s_b1-a_b2-I_0);
figure(1);
freqz(M_l1, 1)
figure(2);
freqz(M_l2,1)
##figure(1);
##plot(t,M_l1)
##hold on
##plot(t,M_l2)
##
##figure(2);
##plot(t,M_b1)
##hold on
##plot(t,M_b2)
##
##figure(1);
##plot(t,s_l1)
##hold on
##plot(t,s_l2)
##hold on
##figure(2);
##plot(t,s_b1)
##hold on
##plot(t,s_b2)
##plot(t,a_l1)
##hold on
##plot(t,a_l2)
##hold on
##plot(t,a_b1)
##hold on
##plot(t,a_b2)
##legend('s_b1', 's_b2', 'position', 'best')
##title('For B_1 and B_2 populations')
%legend('s_l1', 's_l2', 's_b1', 's_b2', 'position', 'best')
##xlim([0,50])
drawnow();

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!