Why am I not getting desired output?
조회 수: 1 (최근 30일)
이전 댓글 표시
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
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
2022년 6월 20일
편집: Torsten
2022년 6월 20일
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
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
2022년 6월 20일
편집: Jan
2022년 6월 20일
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1039120/image.png)
I do not see the problem.
댓글 수: 3
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.
참고 항목
카테고리
Help Center 및 File Exchange에서 Measurements and Spatial Audio에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!