Clear Filters
Clear Filters

Infinite Loop, Values not Updating

1 view (last 30 days)
Irene Zhou
Irene Zhou on 26 Nov 2020
Commented: Irene Zhou on 27 Nov 2020
In the while loop, I am comparing the accuracy of Euler's and Midpoint methods (the methods themselves should be correct). If they differ by a certain amount, I reduce the step size by half and redo the calculation for that time point. Else, I update the indices and move on. Plot(t, x(:,4), t, y(:,4)) should gnerate a graph that portrays an action potential.
Somehow after dif hits the specified threshold, x, y, and dif stop to update. The step size keeps getting reduced and I get stuck in an infinite loop. Does anyone have any insights on what might have happened there? Any thought is appreciated!
dt = 0.08;
Vm = 0;
t(1)=0;
i_inj = 17;
%initial values for n, h, m, V
x = [0.317676914060697 0.596120753508460 0.0529324852572496 0]; %set first value to steady state to 0mV
y = [0.317676914060697 0.596120753508460 0.0529324852572496 0];
i = 1;
k = 1;
while t(k)<11
%Euler's method
x(i+1,:) = x(i,:)+dt.*derivs(t(k),x(i,:),i_inj);
%midpoint method
k1 = dt.*derivs(t(k),y(i,:),i_inj);
k2 = dt.*derivs(t(k)+0.5*dt,y(i,:)+0.5.*k1,i_inj);
y(i+1,:) = k2+y(i,:);
dif = abs(x(i+1,4)-y(i+1,4));
if (dif>=0.1)
dt = dt*0.5;
else
t(k+1) = dt+t(k);
i = i+1;
k = k+1;
end
end
function y = derivs(t,x,I_inj)
g_Kbar = 36;
E_K = -12;
g_Nabar = 120;
E_Na = 115;
Gm = 0.3;
V_rest = 10.613;
V=x(4);
if(V == 10)
alpha_n = 1/10;
else
alpha_n = (10-V)./(100*(exp(1-V/10)-1));
end
beta_n = 0.125*exp(-V/80);
n=x(1);
dndt = alpha_n*(1-n)-beta_n*n;
y(1) = dndt;
alpha_h = 0.07*exp(-V/20);
beta_h = 1/(exp(3-V/10)+1);
h=x(2);
dhdt = alpha_h*(1-h)-beta_h*h;
y(2) = dhdt;
if (V==25)
alpha_m = 1;
else
alpha_m = (25-V)/(10*(exp(2.5-V/10)-1));
end
beta_m = 4*exp(-V/18);
m=x(3);
dmdt = alpha_m*(1-m)-beta_m*m;
y(3) = dmdt;
if t>=1 && t<=1.5
dVdt = I_inj+g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
else
dVdt = g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
end
y(4) = dVdt;
end
  2 Comments
KSSV
KSSV on 26 Nov 2020
What is this function derivs?
Irene Zhou
Irene Zhou on 26 Nov 2020
It's a derivative function that calculates and spits out four variables in a row. This function should be correct as I have tested it in other places.
function y = derivs(t,x,I_inj)
g_Kbar = 36;
E_K = -12;
g_Nabar = 120;
E_Na = 115;
Gm = 0.3;
V_rest = 10.613;
V=x(4);
if(V == 10)
alpha_n = 1/10;
else
alpha_n = (10-V)./(100*(exp(1-V/10)-1));
end
beta_n = 0.125*exp(-V/80);
n=x(1);
dndt = alpha_n*(1-n)-beta_n*n;
y(1) = dndt;
alpha_h = 0.07*exp(-V/20);
beta_h = 1/(exp(3-V/10)+1);
h=x(2);
dhdt = alpha_h*(1-h)-beta_h*h;
y(2) = dhdt;
if (V==25)
alpha_m = 1;
else
alpha_m = (25-V)/(10*(exp(2.5-V/10)-1));
end
beta_m = 4*exp(-V/18);
m=x(3);
dmdt = alpha_m*(1-m)-beta_m*m;
y(3) = dmdt;
if t>=1 && t<=1.5
dVdt = I_inj+g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
else
dVdt = g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
end
y(4) = dVdt;
end

Sign in to comment.

Answers (1)

Andy
Andy on 26 Nov 2020
I don't have real values for n0 etc so I put some random numbers in and the resultant t(k) was generally small, much less than the threshold of 11 that you have set for the while loop. I was just wondering if this is the correct variable for the comparison (should you be using k itself). Can you provide valid examples for your starting values?
  4 Comments
Andy
Andy on 27 Nov 2020
It starts going "wrong" on the 33 run through the loop when t=2.3414, I haven't plotted it but what is the function doing at this point?
Irene Zhou
Irene Zhou on 27 Nov 2020
At the end of my Else statement, I update k = k+1 so t(k) is able to reach 11 in theory.
The purpose of comparing the Euler's method to the midpoint method is to take dynamic step sizes. If the difference is too big, take a smaller step. If the difference is small enough, try a bigger step, which I haven't got to.
It starts to go wrong when dif hits the threshold I set, 0.1 in this case. dt keeps getting reduced by half while nothing else gets updated, which is what I am trying to figure out.

Sign in to comment.

Categories

Find more on App Building in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!