0

Differential equation where M_v is a piecewise function

differential equation where M_v is a piecewise function

Differential equation where M_v is a piecewise function

differential equation where M_v is a piecewise function

Piecewise function

piecewise function

I have coded the differential equation in Octave using the RUnge-Kutta 4 (RK4) method. However, the result is not as desired due to the piecewise function being coded incorrectly in my view. Could anyone help to code piecewise function correctly to solve differential equation using rk4 method?

  %vIRt code of the equations given in the manuscript 

function x = pieceWise(s)
  x = zeros (size (s));
  ind = s > 0;
  x(ind) = s(ind);
endfunction
% Defines variables
global tau_s taur_a beta_r I_ext J_inter J_intra Jr_a  I_0 lp; 

tau_s  = 10; taur_a = 83;  beta_r = 0.0175; Jr_a = 172.97; J_inter= 81; J_intra = 32.4; I_0= 0.29; I_ext = 20;
% step sizes
t0 = 0;
tfinal = 200;
h = 0.05;
n = ceil((tfinal-t0)/h)+1;

%initial conditions
s_r(1) = 0;
s_p(1) = 0.3556;
a_p(1) = 6.151;
a_r(1) = 0;
t(1) =0;
 
% functions handles
S_r = @(t,s_r,s_p,a_r,a_p) -(s_r/tau_s)+ beta_r*pieceWise(I_ext-J_intra*s_r-J_inter*s_p-a_r-I_0);
A_r = @(t,s_r,s_p,a_r,a_p) (-a_r+Jr_a*beta_r*pieceWise(I_ext-J_intra*s_r-J_inter*s_p-a_r-I_0))/taur_a;
S_p = @(t,s_r,s_p,a_p,a_r) -(s_p/tau_s)+ beta_r*pieceWise(I_ext-J_inter*s_r-J_intra*s_p-a_p-I_0);
A_p = @(t,s_r,s_p,a_p,a_r) (-a_p+Jr_a*beta_r*pieceWise((I_ext-J_inter*s_r-J_intra*s_p-a_p-I_0)))/taur_a;


for i = 1:n
  
  %updates time 
  t(i+1) = t(i)+h;
  
  %updates S_r, A_r, S_p and A_p
 k1S_r = S_r(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 k1A_r = A_r(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 k1S_p = S_p(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 k1A_p = A_p(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 
 k2S_r = S_r(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 k2A_r = A_r(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 k2S_p = S_p(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 k2A_p = A_p(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 
 k3S_r = S_r(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 k3A_r = A_r(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 k3S_p = S_p(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 k3A_p = A_p(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 
 k4S_r = S_r(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 k4A_r = A_r(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 k4S_p = S_p(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 k4A_p = A_p(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 
 s_r(i+1) = s_r(i)+h/6*(k1S_r +   2*k2S_r + 2*k3S_r +k4S_r);
 s_p(i+1) = s_p(i)+h/6*(k1S_p +   2*k2S_p + 2*k3S_p +k4S_p);
 a_r(i+1) = a_r(i)+h/6*(k1A_r +   2*k2A_r + 2*k3A_r +k4A_r);
 a_p(i+1) = a_p(i)+h/6*(k1A_p +   2*k2A_p + 2*k3A_p +k4A_p);
end

M_r = beta_r*pieceWise((I_ext-J_intra*s_r-J_inter*s_p-a_r-I_0));
M_p = beta_r*pieceWise((I_ext-J_inter*s_r-J_intra*s_p-a_p-I_0));

%plots
##plot(t,M_p)
##hold on
##plot(t,M_r)
##plot(t,s_r)
##hold on
##plot(t,s_p)
##hold on
##plot(t,a_p)
##hold on
##plot(t,a_r)
##title('a_p and a_r')
##legend('s_r', 's_p', 'position', 'best')
##xlim([0,200])

drawnow;

%print -deps rk4.eps
Adriaan
  • 17,741
  • 7
  • 42
  • 75
Rakesh Kumar
  • 1
  • Please read the descriptions of both [tag:matlab] and [tag:octave]. They are **not** the same, thus please only use both tags when asking about the similarities/differences between the two. – Adriaan May 13 '22 at 09:56
  • Welcome to Stack Overflow! Please take the [tour] and read up on [ask]. Why do you think your piece-wise function is wrong and on which line do you think the problem occurs? – Adriaan May 13 '22 at 09:59
  • When I delete the minus sign from the function handle then still the same plot is showing. And the plot of a_p is constantly increasing with time but the values of right hand side of A_p is always negative that's means a_p should decrease. These are facts indicating that I didn't use piecewise function correctly. – Rakesh Kumar May 13 '22 at 10:13
  • You are aware that in the definition of `S_p` and `A_p` the order of `a_r` and `a_p` is switched? (But that it is not the case in the actual function calls.) – Lutz Lehmann May 13 '22 at 11:16
  • Thank you for your response. I didn't understand your point about how the order of a_r and a_p are switched. Could you explain it in detail? My desired solution is all s_p, s_r, and a_r, a_p should stay at their initial values because they are steady-state values. – Rakesh Kumar May 13 '22 at 16:10
  • The order of inputs you define in your "% function handles" section does not match what you use in the "%updates S_r, A_r, S_p and A_p" section. You switched the order of a_p and a_r when you defined S_p and A_p compared to how you defined the input order in S_r and A_r. If that was not intentional, maybe it's part of the problem. – Nick J May 14 '22 at 13:59
  • Thank you very much Nick and Lehmann for your comments. Now, I am getting the desired results after changing the order of a_r and a_p. – Rakesh Kumar May 15 '22 at 07:26

0 Answers0