Why am I not getting desired output?
1 view (last 30 days)
Show older comments
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
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.
Accepted Answer
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
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;
More Answers (1)
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
Jan
on 20 Jun 2022
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.
See Also
Categories
Find more on Measurements and Spatial Audio in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!