Why am I not getting desired output?

3 次查看（过去 30 天）
Rakesh
Rakesh 2022-6-20
评论： Torsten 2022-6-20
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 个评论
Jan
Jan 2022-6-20
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.

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2022-6-20
编辑：Torsten 2022-6-20
Ran in:
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 个评论
Rakesh
Rakesh 2022-6-20
Thank you for your support!
Could you let know where did I do mistake?
Torsten
Torsten 2022-6-20
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;

请先登录，再进行评论。

更多回答（1 个）

Jan
Jan 2022-6-20
编辑：Jan 2022-6-20
Ran in:
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 个评论
显示 1更早的评论
Jan
Jan 2022-6-20
Please copy and paste the code from my answer to your Octave and run it.
I'm convinced, that you use a different code - e.g. your diagram ends and 200.
Rakesh
Rakesh 2022-6-20
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();

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Octave 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by