Why am I not getting desired output?
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commento
Jan
il 20 Giu 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.
Risposta accettata
Torsten
il 20 Giu 2022
Modificato: Torsten
il 20 Giu 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 Commenti
Torsten
il 20 Giu 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;
Più risposte (1)
Jan
il 20 Giu 2022
Modificato: Jan
il 20 Giu 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1039120/image.png)
I do not see the problem.
3 Commenti
Jan
il 20 Giu 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.
Vedere anche
Categorie
Scopri di più su Measurements and Spatial Audio in Help Center e File Exchange
Prodotti
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!